1
$\begingroup$

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?

$\endgroup$

1 Answer 1

0
$\begingroup$

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)
  
}
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.