R package to solve Gaussian MLE under conditional independence constraints

Question

Is there any R package or function to solve Gaussian MLE under conditional independence constraints?

Suppose we have $y_i\overset{i.i.d}{\sim}\mathcal{N}(0,\Sigma_{p\times p})$, $i = 1,2,\ldots,n$. We know that $(\Sigma^{-1})_{ij} = 0$, for some $i,j\in\{1,2,\ldots,p\}$, i.e. $X_i$ and $X_j$ are conditionally independent given the rest of the variables. We would like to find MLE of $\Sigma$ under the conditional independence constraints.

I also tried to implement it according to section 5.1 of the above paper, but I couldn't do it successfully.

I wonder if there is any R implementation to find Gaussian MLE under conditional independence constraints?

score 0 · Answer 1 · answered Oct 08 '21 at 15:17

I finished the implementation myself. Here is the code.

gaussian_MLE = function(S, omega_structure, initial_omega = initial_omega, maxit = 2e9, tol = 1e-15){
  omegas = list()
  
  
  # S is sample covariance matrix; omega_structure is the adjacency matrix. 
  
  temp = which(omega_structure == 1, arr.ind = TRUE)
  temp = temp[temp[,1]>=temp[,2],]
  
  r = 1
  converged = FALSE
  
  omega_old = initial_omega
  omega_old[!omega_structure] = 0
  omega_current = omega_old
  while(!converged){
   
    
    omega_old = omega_current
    
    for(i in 1:nrow(temp)){

      S_current = solve(omega_current)
      
      u = temp[i,1]
      v = temp[i,2]
      
      sqrt = sqrt(S_current[u,u] * S_current[v,v])
      
      rho_bar = S[u, v]/sqrt
      # print(rho_bar)
      
      rho = S_current[u, v]/sqrt
      
      if(rho==1){
        if(rho_bar > 0){
          s = (1-rho_bar)/(2*rho_bar)/sqrt
        }else{
          stop("here1")
        }
      }else if(rho == -1){
        if(rho_bar < 0){
          s = (1+rho_bar)/(2*rho_bar)/sqrt
        }else{
          stop("here2")
        }
      }else{
        one_minus_rho_square = 1-rho^2
        a = rho_bar*one_minus_rho_square
        b = -(one_minus_rho_square+2*rho*rho_bar)
        c = rho-rho_bar
        
        if(rho_bar==0){
          s = rho/one_minus_rho_square/sqrt
        }else if(rho_bar>0){
          s = min((-b+sqrt(b^2-4*a*c))/(2*a)/sqrt,(-b-sqrt(b^2-4*a*c))/(2*a)/sqrt)
        }else{
          s = max((-b+sqrt(b^2-4*a*c))/(2*a)/sqrt,(-b-sqrt(b^2-4*a*c))/(2*a)/sqrt)
        }
        
      }
      
      if(u==v){
        omega_current[u,v] =  omega_current[u,v] +  s
      }else{
        omega_current[u,v] =  omega_current[u,v] +  s
        omega_current[v,u] =  omega_current[v,u] +  s
      }
      

      
    }
    
    if(max(abs(omega_current - omega_old))<tol){
      converged = TRUE
    }else{
      r = r+1
      if((r %% 1000000)==0){
        print(r)
      }
    }
    
    if(r>maxit){
      stop("Algorithm did not converge.")
      break;
    }
    
    # print(r)
    # print(omega_old)

    # omegas[[length(omegas)+1]] = omega_old
    # omegas <<- omegas
    
  }
  
  
  if(max(abs(solve(omega_current)*omega_structure - S*omega_structure))>1e-13){
    stop("here3")
  }
  eigen(omega_current)
  
  
  
  rrr <<- r
  
  return(omega_current)
  
}

R package to solve Gaussian MLE under conditional independence constraints

1 Answers1