Move all fitfunction args to expectation
[openmx:openmx.git] / models / nightly / BootstrapParallelNightly.R
1 #
2 #   Copyright 2007-2010 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 require(OpenMx)
17 require(snowfall)
18 sfInit(parallel = TRUE, cpus = 8)
19 sfLibrary(OpenMx)
20
21 set.seed(10)
22
23 # parameters for the simulation: lambda = factor loadings,
24 # specifics = specific variances
25 lambda <- matrix(c(.8, .5, .7, 0), 4, 1)
26 nObs <- 500
27 nReps <- 1000
28 nVar <- nrow(lambda)
29 specifics <- diag(nVar)
30 chl <- chol(lambda %*% t(lambda) + specifics)
31
32 # indices for parameters and hessian estimate in results
33 pStrt <- 3
34 pEnd <- pStrt + 2*nVar - 1
35 hStrt <- pEnd + 1
36 hEnd <- hStrt + 2*nVar - 1
37
38 # dimension names for OpenMx
39 dn <- list()
40 dn[[1]] <- paste("Var", 1:4, sep="")
41 dn[[2]] <- dn[[1]]
42
43 # function to get a covariance matrix
44 randomCov <- function(nObs, nVar, chl, dn) {
45   x <- matrix(rnorm(nObs*nVar), nObs, nVar)
46   x <- x %*% chl
47   thisCov <- cov(x)
48   dimnames(thisCov) <- dn
49   return(thisCov)  
50 }
51
52 createNewModel <- function(model, modelname) {
53         model@data@observed <- randomCov(nObs, nVar, chl, dn)
54         model@name <- modelname
55         return(model)
56 }
57
58 getStats <- function(model) {
59         retval <- c(model@output$status[[1]],
60                 max(abs(model@output$gradient)),
61                 model@output$estimate,
62                 sqrt(diag(solve(model@output$hessian))))
63         return(retval)
64 }
65
66
67 # initialize obsCov for MxModel
68 obsCov <- randomCov(nObs, nVar, chl, dn)
69
70 # results matrix: get results for each simulation
71 results <- matrix(0, nReps, hEnd)
72 dnr <- c("inform", "maxAbsG", paste("lambda", 1:nVar, sep=""),
73          paste("specifics", 1:nVar, sep=""),
74          paste("hessLambda", 1:nVar, sep=""),
75          paste("hessSpecifics", 1:nVar, sep=""))
76 dimnames(results)[[2]] <- dnr
77
78 # instantiate MxModel
79 template <- mxModel(name="stErrSim",
80                        mxMatrix(name="lambda", type="Full", nrow=4, ncol=1,
81                                 free=TRUE, values=c(.8, .5, .7, 0)),
82                        mxMatrix(name="specifics", type="Diag", nrow=4,
83                                 free=TRUE, values=rep(1, 4)),
84                        mxAlgebra(lambda %*% t(lambda) + specifics,
85                                  name="preCov", dimnames=dn),
86                        mxData(observed=obsCov, type="cov", numObs=nObs),
87                        mxFitFunctionML(),mxExpectationNormal(covariance='preCov'),
88                        independent = TRUE)
89
90 topModel <- mxModel(name = 'container')
91
92 submodels <- lapply(1:nReps, function(x) {
93   createNewModel(template, paste('stErrSim', x, sep=''))
94 })
95
96 names(submodels) <- imxExtractNames(submodels)
97 topModel@submodels <- submodels
98
99 modelResults <- mxRun(topModel, silent=TRUE, suppressWarnings=TRUE)
100
101 results <- t(omxSapply(modelResults@submodels, getStats))
102
103 sfStop()
104
105 # get rid of bad covergence results
106 results2 <- data.frame(results[which(results[,1] <= 1),])
107
108 # summarize the results
109 means <- colMeans(results2)
110 stdevs <- sapply(results2, sd)
111 sumResults <- data.frame(matrix(dnr[pStrt:pEnd], 2*nVar, 1,
112                                 dimnames=list(NULL, "Parameter")))
113 sumResults$mean <- means[pStrt:pEnd]
114 sumResults$obsStDev <- stdevs[pStrt:pEnd]
115 sumResults$meanHessEst <- means[hStrt:hEnd]
116 sumResults$sqrt2meanHessEst <- sqrt(2) * sumResults$meanHessEst
117
118 omxCheckCloseEnough(sumResults[1,'mean'], 0.8, 0.1)
119 omxCheckCloseEnough(sumResults[2,'mean'], 0.5, 0.1)
120 omxCheckCloseEnough(sumResults[3,'mean'], 0.7, 0.1)
121 omxCheckCloseEnough(sumResults[4,'mean'], 0.0, 0.1)
122 omxCheckCloseEnough(sumResults[5,'mean'], 1.00, 0.1)
123 omxCheckCloseEnough(sumResults[6,'mean'], 1.00, 0.1)
124 omxCheckCloseEnough(sumResults[7,'mean'], 1.00, 0.1)
125 omxCheckCloseEnough(sumResults[8,'mean'], 1.00, 0.1)
126
127 detach("package:snowfall")
128 unloadNamespace("snowfall")
129 detach("package:snow")