EM algorithm for incomplete normal data


Performs maximum-likelihood estimation on the matrix of incomplete data using the EM algorithm. Can also be used to find a posterior mode under a normal-inverted Wishart prior supplied by the user.


em.norm(s, start, showits=TRUE, maxits=1000, criterion=0.0001, prior)



summary list of an incomplete normal data matrix produced by the function prelim.norm.


optional starting value of the parameter. This is a parameter vector in packed storage, such as one created by the function makeparam.norm. If no starting value is supplied, em.norm chooses its own starting value.


if TRUE, reports the iterations of EM so the user can monitor the progress of the algorithm.


maximum number of iterations performed. The algorithm will stop if the parameter still has not converged after this many iterations.


convergence criterion. The algorithm stops when the maximum relative difference in all of the estimated means, variances, or covariances from one iteration to the next is less than or equal to this value.


optional prior distribution. This is a list containing the hyperparameters of a normal-inverted Wishart distribution. In order, the elements of the list are: tau (a scalar), m (a scalar), mu0 (a vector of length ncol(x)), and lambdainv (a matrix of dimension c(ncol(x),ncol(x))). The elements of mu0 ans lambdainv apply to the data after transformation, i.e. after the columns have been centered and scaled to have mean zero and variance one. If no prior is supplied, the default is a uniform prior, which results in maximum-likelihood estimation.


The default starting value takes all means on the transformed scale to be equal to zero, and covariance matrix on the transformed scale equal to the identity. All important computations are carried out in double precision, using the sweep operator.


a vector representing the maximum-likelihood estimates of the normal parameters. This vector contains means, variances, and covariances on the transformed scale in packed storage. The parameter can be transformed back to the original scale and put into a more understandable format by the function getparam.norm.


See Section 5.3 of Schafer (1994).

See Also

prelim.norm, getparam.norm, and makeparam.norm.


s <- prelim.norm(mdata)   #do preliminary manipulations
thetahat <- em.norm(s)   #compute mle
getparam.norm(s,thetahat,corr=TRUE)$r #look at estimated correlations

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker. Vote for new features on Trello.

comments powered by Disqus