mxComputeEM {OpenMx}

Fit a model using DLR's (1977) Expectation-Maximization (EM) algorithm

Description

The EM algorithm constitutes the following steps: Start with an initial parameter vector. Predict the missing data to form a completed data model. Optimize the completed data model to obtain a new parameter vector. Repeat these steps until convergence criteria are met.

Usage

mxComputeEM(
  expectation = NULL,
  predict = NA_character_,
  mstep,
  observedFit = "fitfunction",
  ...,
  maxIter = 500L,
  tolerance = 1e-09,
  verbose = 0L,
  freeSet = NA_character_,
  accel = "varadhan2008",
  information = NA_character_,
  infoArgs = list(),
  estep = NULL
)

Arguments

expectation

a vector of expectation names

predict

what to predict from the observed data

mstep

a compute plan to optimize the completed data model

observedFit

the name of the observed data fit function (defaults to "fitfunction")

...

Not used. Forces remaining arguments to be specified by name.

maxIter

maximum number of iterations

tolerance

optimization is considered converged when the maximum relative change in fit is less than tolerance

verbose

integer. Level of run-time diagnostic output. Set to zero to disable

freeSet

names of matrices containing free variables

accel

name of acceleration method ("varadhan2008" or "ramsay1975")

information

name of information matrix approximation method

infoArgs

arguments to control the information matrix method

estep

a compute plan to perform the expectation step

Details

The arguments to this function have evolved. The old style mxComputeEM(e,p,mstep=m) is equivalent to the new style mxComputeEM(estep=mxComputeOnce(e,p), mstep=m). This change allows the API to more closely match the literature on the E-M method. You might use mxAlgebra(..., recompute='onDemand') to contain the results of the E-step and then cause this algebra to be recomputed using mxComputeOnce.

This compute plan does not work with any and all expectations. It requires a special kind of expectation that can predict its missing data to create a completed data model.

The EM algorithm does not produce a parameter covariance matrix for standard errors. The Oakes (1999) direct method and S-EM, an implementation of Meng & Rubin (1991), are included.

Ramsay (1975) was recommended in Bock, Gibbons, & Muraki (1988).

References

Bock, R. D., Gibbons, R., & Muraki, E. (1988). Full-information item factor analysis. Applied Psychological Measurement, 6(4), 431-444.

Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society. Series B (Methodological), 1-38.

Meng, X.-L. & Rubin, D. B. (1991). Using EM to obtain asymptotic variance-covariance matrices: The SEM algorithm. Journal of the American Statistical Association, 86 (416), 899-909.

Oakes, D. (1999). Direct calculation of the information matrix via the EM algorithm. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 61(2), 479-482.

Ramsay, J. O. (1975). Solving implicit equations in psychometric data analysis. Psychometrika, 40 (3), 337-360.

Varadhan, R. & Roland, C. (2008). Simple and globally convergent methods for accelerating the convergence of any EM algorithm. Scandinavian Journal of Statistics, 35, 335-353.

See Also

MxAlgebra, mxComputeOnce

Examples

library(OpenMx)
set.seed(190127)

N <- 200
x <- matrix(c(rnorm(N/2,0,1),
              rnorm(N/2,3,1)),ncol=1,dimnames=list(NULL,"x"))
data4mx <- mxData(observed=x,type="raw")

class1 <- mxModel("Class1",
	mxMatrix(type="Full",nrow=1,ncol=1,free=TRUE,values=0,name="Mu"),
	mxMatrix(type="Full",nrow=1,ncol=1,free=TRUE,values=4,name="Sigma"),
	mxExpectationNormal(covariance="Sigma",means="Mu",dimnames="x"),
	mxFitFunctionML(vector=TRUE))

class2 <- mxRename(class1, "Class2")

mm <- mxModel(
	"Mixture", data4mx, class1, class2,
	mxAlgebra((1-Posteriors) * Class1.fitfunction, name="PL1"),
	mxAlgebra(Posteriors * Class2.fitfunction, name="PL2"),
	mxAlgebra(PL1 + PL2, name="PL"),
	mxAlgebra(PL2 / PL,  recompute='onDemand',
	          initial=matrix(runif(N,.4,.6), nrow=N, ncol = 1), name="Posteriors"),
	mxAlgebra(-2*sum(log(PL)), name="FF"),
	mxFitFunctionAlgebra(algebra="FF"),
	mxComputeEM(
	  estep=mxComputeOnce("Mixture.Posteriors"),
	  mstep=mxComputeGradientDescent(fitfunction="Mixture.fitfunction")))

mm <- mxOption(mm, "Max minutes", 1/20)  # remove this line to find optimum
mmfit <- mxRun(mm)
summary(mmfit)