copied scripts from passing
[openmx:openmx.git] / demo / BivariateSaturated_MatrixRawCholesky.R
1 require(OpenMx)
2
3 #Simulate Data
4 require(MASS)
5 set.seed(200)
6 rs=.5
7 xy <- mvrnorm (1000, c(0,0), matrix(c(1,rs,rs,1),2,2))
8 testData <- xy
9 selVars <- c('X','Y')
10 dimnames(testData) <- list(NULL, selVars)
11 summary(testData)
12 cov(testData)
13
14 #examples 6: Saturated Model with RawData and Matrix-Style Input
15 bivSatModel6 <- mxModel("bivSat6",
16     mxMatrix(
17         type="Full", 
18         nrow=2, 
19         ncol=2, 
20         free=c(T,T,F,T), 
21         values=c(1,.2,0,1), 
22         name="Chol"
23     ),
24     mxAlgebra(
25         Chol %*% t(Chol), 
26         name="expCov", 
27         dimnames=list(selVars,selVars)
28     ),
29     mxMatrix(
30         type="Full", 
31         nrow=1, 
32         ncol=2, 
33         free=T, 
34         values=c(0,0), 
35         dimnames=list(NULL,selVars), 
36         name="expMean"
37     ),
38     mxData(
39         observed=testData, 
40         type="raw", 
41     ),
42     mxFIMLObjective(
43         covariance="expCov",
44         means="expMean"
45     )
46     )
47 bivSatFit6 <- mxRun(bivSatModel6)
48 EM6 <- mxEvaluate(expMean, bivSatFit6)
49 EC6 <- mxEvaluate(expCov, bivSatFit6)
50 LL6 <- mxEvaluate(objective,bivSatFit6)
51
52
53 #Mx answers hard-coded
54 #example Mx..2: Saturated Model with Raw Data
55 Mx.EM2 <- matrix(c(0.03211188, -0.004889211),1,2)
56 Mx.EC2 <- matrix(c(1.0092891, 0.4813504, 0.4813504, 0.9935366),2,2)
57 Mx.LL2 <- 5415.772
58
59
60 #Compare OpenMx results to Mx results (LL: likelihood; EC: expected covariance, EM: expected means)
61 #6:RawMat Cholesky
62 omxCheckCloseEnough(LL6,Mx.LL2,.001)
63 omxCheckCloseEnough(EC6,Mx.EC2,.001)
64 omxCheckCloseEnough(EM6,Mx.EM2,.001)
65