Integrating PPML analytical solution in to the trunk.
[openmx:openmx.git] / models / passing / PPML_100days.R
1 # One hundred days model\r
2 # Each manifest represents a single variable as measured on a different day, days 1-100\r
3 # Three predicting latents. Over the trial, one represents a constant part,\r
4 # one represents a linear part, and one an exponential part\r
5 \r
6 require(OpenMx)\r
7 require(MASS)\r
8 \r
9 # Random covariance matrix\r
10 #dataTest <- matrix(rnorm(100*100), 100, 100)\r
11 #dataTest <- dataTest %*% t(dataTest)\r
12 \r
13 manifests <- unlist(lapply(1:100, function(n) { paste("Day",n, sep="") } ))\r
14 latents <- c('Const', 'Lin', 'Exp')\r
15 \r
16 amtT <- 100\r
17 amtStep <- amtT / 99\r
18 ConstLoadings <- rep(1,100)\r
19 LinLoadings <- unlist(lapply(1:100, function (n) { n*amtStep/(amtT-1) } ))\r
20 EXPFACTOR <- 0.22\r
21 ExpLoadings <- unlist(lapply(1:100, function (n) { -exp(-EXPFACTOR*n*amtStep) } ))\r
22 \r
23 lambda <- cbind(ConstLoadings, LinLoadings, ExpLoadings)\r
24 dataLatents <- matrix(rnorm(3*3), 3, 3)\r
25 dataTest <- lambda %*% dataLatents %*% t(dataLatents) %*% t(lambda) + diag(rep(1,100))\r
26 dataTest <- (dataTest + t(dataTest)) / 2\r
27 \r
28 dataRaw <- mvrnorm(n=5000, rep(0,100), dataTest)\r
29 colnames(dataRaw) <- manifests\r
30 \r
31 browser()\r
32 \r
33 colnames(dataTest) <- manifests\r
34 rownames(dataTest) <- manifests\r
35 \r
36 factorModel <- mxModel("One Hundred Days Model",\r
37       type="RAM",\r
38           # Vars\r
39       manifestVars = manifests,\r
40       latentVars = latents,\r
41       # A\r
42           mxPath(from='Const',  to=manifests,value=ConstLoadings,       free=FALSE),\r
43           mxPath(from='Lin',    to=manifests,value=LinLoadings,         free=FALSE),\r
44           mxPath(from='Exp',    to=manifests,value=ExpLoadings,         free=FALSE),\r
45       # S\r
46           mxPath(from=manifests, arrows=2,value=rep(1.0, 100), labels=rep("Res", 100)),\r
47           #mxPath(from=latents, arrows=2,values=1.0),\r
48           mxPath(from="Const", to=latents, arrows=2,values=1.0),\r
49           mxPath(from="Lin", to=latents, arrows=2,values=1.0),\r
50           mxPath(from="Exp", to=latents, arrows=2,values=1.0),\r
51           \r
52           # Means\r
53           mxPath(from="one", to=latents, values=0, free=TRUE),\r
54           \r
55           # Data\r
56       #mxData(dataTest, type="cov", numObs=100))\r
57           mxData(dataRaw, type="raw", numObs=5000)\r
58         )\r
59 \r
60 imxTestPPML(factorModel)\r