old files deleted, new files added
[openmx:openmx.git] / demo / BivariateSaturated_MatrixRawCholesky.R
1 # -----------------------------------------------------------------------
2 # Program: BivariateSaturated_MatrixRawCholesky.R  
3 #  Author: Hermine Maes
4 #    Date: 08 01 2009 
5 #
6 # Bivariate Saturated model to estimate means and (co)variances 
7 # using Cholesky Decomposition
8 # Matrix style model input - Raw data input
9 #
10 # Revision History
11 #   Hermine Maes -- 10 08 2009 updated & reformatted
12 # -----------------------------------------------------------------------
13
14 require(OpenMx)
15
16 #Simulate Data
17 # -----------------------------------------------------------------------
18 require(MASS)
19 set.seed(200)
20 rs=.5
21 xy <- mvrnorm (1000, c(0,0), matrix(c(1,rs,rs,1),2,2))
22 testData <- xy
23 selVars <- c("X","Y")
24 dimnames(testData) <- list(NULL, selVars)
25 summary(testData)
26 cov(testData)
27
28 #example 6: Saturated Model with Raw Data and Matrix-Style Input
29 # -----------------------------------------------------------------------
30 bivSatModel6 <- mxModel("bivSat6",
31     mxMatrix(
32                 type="Full", 
33                 nrow=2, 
34                 ncol=2, 
35                 free=c(TRUE,TRUE,FALSE,TRUE), 
36                 values=c(1,.2,0,1), 
37                 name="Chol"
38         ),
39         mxMatrix(
40                 type="Full", 
41                 nrow=1, 
42                 ncol=2, 
43                 free=TRUE, 
44                 values=c(0,0), 
45                 name="expMean"
46         ),
47     mxAlgebra(
48                 expression=Chol %*% t(Chol), 
49                 name="expCov"
50         ),
51         mxData(
52                 observed=testData, 
53                 type="raw"
54         ),
55         mxFIMLObjective(
56                 covariance="expCov",
57                 means="expMean",
58                 dimnames=selVars
59         )
60 )
61
62 bivSatFit6 <- mxRun(bivSatModel6)
63 EM <- mxEval(expMean, bivSatFit6)
64 EC <- mxEval(expCov, bivSatFit6)
65 LL <- mxEval(objective,bivSatFit6)
66
67 #Mx answers hard-coded
68 # -----------------------------------------------------------------------
69 #example Mx..2: Saturated Model with Raw Data
70 Mx.EM <- matrix(c(0.03211188, -0.004889211),1,2)
71 Mx.EC <- matrix(c(1.0092891, 0.4813504, 0.4813504, 0.9935366),2,2)
72 Mx.LL <- 5415.772
73
74 # Compare OpenMx results to Mx results 
75 # -----------------------------------------------------------------------
76 # (LL: likelihood; EC: expected covariance, EM: expected means)
77 #6: RawMat Cholesky
78 omxCheckCloseEnough(LL,Mx.LL,.001)
79 omxCheckCloseEnough(EC,Mx.EC,.001)
80 omxCheckCloseEnough(EM,Mx.EM,.001)