title: Running error within function using solve() in R?
body: <p>With this program below, I will get the error:</p>
  
  <pre><code>solve.default(Sigma0[cs.idx, cs.idx]) : 'a' is 0-diml
  </code></pre>
  
  <p>But, when I check the <code>em()</code> function step by step, I mean, sentence by sentence without function, there is no error within <code>solve()</code>. So I am confused and desperate for help, Thank you!</p>
  
  <pre><code>###----------------------------------------------------------------
  ### Maximal Likelihood estimation of mean and covariance 
  ###  for multivariate normal distribution by EM algorithm,
  ###             for demonstration purposes only
  ###----------------------------------------------------------------
  
  em&lt;-function(xdata,mu0,Sigma0){
    n&lt;-nrow(xdata)
    p&lt;-ncol(xdata)
  
     err&lt;-function(mu0,Sigma0,mu1,Sigma1){
      th0&lt;-c(mu0,as.vector(Sigma0))
      th1&lt;-c(mu1,as.vector(Sigma1))
      sqrt(sum((th0-th1)*(th0-th1)))
     }
  
    mu1&lt;-mu0+1
    Sigma1&lt;-Sigma0+1
    while(err(mu0,Sigma0,mu1,Sigma1)&gt;1e-6){
       mu1&lt;-mu0
       Sigma1&lt;-Sigma0
  
       zdata&lt;-xdata
       Ai&lt;-matrix(0,p,p)
      for(i in 1:n){
       if(any(is.na(xdata[i,]))){
        zi&lt;-xdata[i,]
        na.idx&lt;-(1:p)[is.na(zi)]
        cs.idx&lt;-(1:p)[-na.idx]
        Sigma012&lt;-Sigma0[na.idx,cs.idx,drop=FALSE]
        Sigma022.iv&lt;-solve(Sigma0[cs.idx,cs.idx])
        zdata[i,na.idx]&lt;-mu0[na.idx]+(Sigma012%*%Sigma022.iv)%*%(zi[cs.idx]-mu0[cs.idx])
        Ai[na.idx,na.idx]&lt;-Ai[na.idx,na.idx]+Sigma0[na.idx,na.idx]-Sigma012%*%Sigma022.iv%*%t(Sigma012)
         }
       }
      mu0&lt;-colMeans(zdata)
      Sigma0&lt;-(n-1)*cov(zdata)/n+Ai/n
     }
   return(list(mu=mu0,Sigma=Sigma0))
  }
  
  ##A simulation example
  library(MASS)
  set.seed(1200)
  p=3
  mu&lt;-c(1,0,-1)
  n&lt;-1000
  Sig &lt;- matrix(c(1, .7, .6, .7, 1, .4, .6, .4, 1), nrow = 3)
  triv&lt;-mvrnorm(n,mu,Sig)
  misp&lt;-0.2   #MCAR probability
  misidx&lt;-matrix(rbinom(3*n,1,misp)==1,nrow=n)
  triv[misidx]&lt;-NA 
  
  #exclude the cases whose entire elements were missed
  er&lt;-which(apply(apply(triv,1,is.na),2,sum)==p)   
  if(length(er)&gt;=1) triv&lt;-triv[-er,]
  
  #initial values
  mu0&lt;-rep(0,p)
  Sigma0&lt;-diag(p)
  system.time(rlt&lt;-em(triv,mu0,Sigma0))
  
  #a better initial values
  mu0&lt;-apply(triv,2,mean,na.rm=TRUE)
  nas&lt;-is.na(triv)
  na.num&lt;-apply(nas,2,sum)
  zdata&lt;-triv
  zdata[nas]&lt;-rep(mu0,na.num)
  Sigma0&lt;-cov(zdata)
  
  system.time(rlt&lt;-em(triv,mu0,Sigma0))
  </code></pre>