Updated copyright to 2013 for R/ demo/ models/passing and src/ folders, and also...
[openmx:openmx.git] / demo / BivariateSaturated_MatrixRawCholesky.R
1 #
2 #   Copyright 2007-2013 The OpenMx Project
3 #
4 #   Licensed under the Apache License, Version 2.0 (the "License");
5 #   you may not use this file except in compliance with the License.
6 #   You may obtain a copy of the License at
7
8 #        http://www.apache.org/licenses/LICENSE-2.0
9
10 #   Unless required by applicable law or agreed to in writing, software
11 #   distributed under the License is distributed on an "AS IS" BASIS,
12 #   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 #   See the License for the specific language governing permissions and
14 #   limitations under the License.
15
16
17 # -----------------------------------------------------------------------------
18 # Program: BivariateSaturated_MatrixRawCholesky.R  
19 # Author: Hermine Maes
20 # Date: 2009.08.01 
21 #
22 # ModelType: Saturated
23 # DataType: Continuous
24 # Field: None
25 #
26 # Purpose: 
27 #      Bivariate Saturated model to estimate means and (co)variances 
28 #      using Cholesky Decomposition
29 #      Matrix style model input - Raw data input
30 #
31 # RevisionHistory:
32 #      Hermine Maes -- 2009.10.08 updated & reformatted
33 #      Ross Gore -- 2011.06.15 added Model, Data & Field metadata     
34 # -----------------------------------------------------------------------------
35
36 require(OpenMx)
37 require(MASS)
38 # Load Libraries
39 # -----------------------------------------------------------------------------
40
41 set.seed(200)
42 rs=.5
43 xy <- mvrnorm (1000, c(0,0), matrix(c(1,rs,rs,1),2,2))
44 testData <- xy
45 selVars <- c("X","Y")
46 dimnames(testData) <- list(NULL, selVars)
47 summary(testData)
48 cov(testData)
49 # Simulate Data
50 # -----------------------------------------------------------------------------
51
52 bivSatModel6 <- mxModel("bivSat6",
53     mxMatrix(
54                 type="Full", 
55                 nrow=2, 
56                 ncol=2, 
57                 free=c(TRUE,TRUE,FALSE,TRUE), 
58                 values=c(1,.2,0,1), 
59                 name="Chol"
60         ),
61         mxMatrix(
62                 type="Full", 
63                 nrow=1, 
64                 ncol=2, 
65                 free=TRUE, 
66                 values=c(0,0), 
67                 name="expMean"
68         ),
69     mxAlgebra(
70                 expression=Chol %*% t(Chol), 
71                 name="expCov"
72         ),
73         mxData(
74                 observed=testData, 
75                 type="raw"
76         ),
77         mxFitFunctionML(),mxExpectationNormal(
78                 covariance="expCov",
79                 means="expMean",
80                 dimnames=selVars
81         )
82 )
83
84 bivSatFit6 <- mxRun(bivSatModel6)
85 EM <- mxEval(expMean, bivSatFit6)
86 EC <- mxEval(expCov, bivSatFit6)
87 LL <- mxEval(objective,bivSatFit6)
88 # example 6: Saturated Model with Raw Data and Matrix-Style Input
89 # -----------------------------------------------------------------------------
90
91 Mx.EM <- matrix(c(0.03211188, -0.004889211),1,2)
92 Mx.EC <- matrix(c(1.0092891, 0.4813504, 0.4813504, 0.9935366),2,2)
93 Mx.LL <- 5415.772
94 # Mx answers hard-coded
95 # -------------------------------------
96 # example Mx..2: Saturated Model with Raw Data
97 # -----------------------------------------------------------------------------
98
99 omxCheckCloseEnough(LL,Mx.LL,.001)
100 omxCheckCloseEnough(EC,Mx.EC,.001)
101 omxCheckCloseEnough(EM,Mx.EM,.001)
102 # 6: RawMat Cholesky
103 # -------------------------------------
104 # Compare OpenMx results to Mx results 
105 # (LL: likelihood; EC: expected covariance, EM: expected means)
106 # -----------------------------------------------------------------------------