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