# load or create data

# sample size
n <- 200
set.seed(10)
fac <- matrix(rnorm(n), n, 1)
load <- matrix(.8, 1, 4)
e <- matrix(rnorm(4*n, 0, .6), n, 4)
sim1 <- fac %*% load + e
#obs <- round(data.frame(sim1), 2)
obs <- data.frame(matrix(0, n, 4))
obs[sim1>0] <- 1

obs <- mxFactor(data.frame(obs), c(0,1))

names(obs) <- paste("x", 1:4, sep="")

# run the model
# EVERY PARAMETER MUST HAVE A LABEL!!!

model <- mxModel("Alg",
	mxData(obs, "raw"),
	mxMatrix("Full", 4, 1, TRUE, 1, 
		paste("l", 1:4, sep=""), name="A"),
	mxMatrix("Diag", 4, 4, FALSE, 0,
		paste("E[", 1:4, ",1]", sep=""), name="U"),
	mxAlgebra(A %*% t(A) + U, name="cov"),
	mxMatrix("Zero", 1, 4, name="M"),
	mxMatrix("Full", 1, 4, TRUE, 0,
		paste("t", 1:4, sep=""), 
		dimnames=list(NA, names(obs)),
		name="T"),
	mxMatrix("Unit", 4, 1, name="J"),
	mxAlgebra(J-A%^%2, name="E"),
	mxFIMLObjective("cov", "M", names(obs), "T")
	)
	
first <- mxRun(model)	

first@output$calculatedHessian