Latent distribution parameters, gradient and information
[openmx:openmx.git] / models / passing / ifa-latent-2tier.R
1 # echo 0 > /proc/self/coredump_filter  # normally 023
2 # R --vanilla --no-save -f models/failing/bock-aitkin-1981.R
3 # R -d gdb --vanilla --no-save -f models/failing/bock-aitkin-1981.R
4
5 #options(error = browser)
6 require(OpenMx)
7 require(rpf)
8 library(mvtnorm)
9
10 set.seed(1)
11
12 numItems <- 16
13 numPeople <- 500
14
15 items <- list()
16 correct <- list()
17 for (ix in 1:numItems) {
18   if (ix <= 5 || ix > numItems-5) {
19     items[[ix]] <- rpf.grm(factors=2)
20   } else {
21     items[[ix]] <- rpf.grm(factors=3)
22   }
23         correct[[ix]] <- rpf.rparam(items[[ix]])
24 }
25 maxParam <- max(vapply(items, function(i) rpf.numParam(i), 0))
26 correct.mat <- sapply(correct, function(p) c(p, rep(NA, maxParam-length(p))))
27
28 true.mean <- c(.2, -.1, .3, -.2, -.5, .15)
29 true.cov <- diag(6)
30 diag(true.cov) <- c(.5, 1.2, 1.5, .9, 1.1, 1.3)
31 true.cov[1,2] <- -.2
32 true.cov[2,1] <- -.2
33 true.pt <- c(true.mean, diag(true.cov)[-2], true.cov[2,1])
34
35 design <- matrix(c(rep(1L,numItems-5), rep(NA,5),
36                    rep(NA,5), rep(2L, numItems-5),
37                    as.integer(kronecker(3:6,rep(1,4)))), byrow=TRUE, ncol=numItems)
38
39 data <- rpf.sample(numPeople, items, correct.mat, design=design, mean=true.mean, cov=true.cov)
40
41 ip.mat <- mxMatrix(name="ItemParam", nrow=maxParam, ncol=numItems,
42                    values=correct.mat, free=FALSE)
43
44 cache <- list(c(-154.47824, 68.38876, -141.9051, 168.81419, -163.4029, 70.22361,  
45                 34.2317, 30.47046, -67.53613, -20.37652, -17.7683, 3.54871, 64.94546 ),
46               c(-513.74756, 265.22542, -300.76141, -252.1387, 515.74301, 134.22559,  
47                 -163.65775, 14.80394, -106.28184, -28.84007, -184.27976, -17.49414,  180.68775),
48               c(-546.05265, 330.52021, -382.72443, -58.10493, 231.0593, 125.1684,  
49                 -153.4536, -40.19255, -83.96658, 10.302, -10.73235, -0.7226,  246.54888),
50               c(107.54376, -315.77402, -174.51499, -71.11526, 238.98663, -303.64933,  
51                 -20.35021, -80.25071, -69.56831, -64.78418, -12.00931, -81.26354,  208.5762),
52               c(-323.44541, 225.38161, 52.08564, -109.818, -92.97802, 263.88848,  
53                 -57.49684, -26.74753, -36.9511, -26.04965, 9.54675, -65.79383,  229.31144),
54               c(7.90011, 442.91977, -188.27667, -51.21122, 440.24953, 147.90806,  
55                 55.09937, -52.71043, -8.95109, -51.11931, -117.02963, -11.38825,  36.38522),
56               c(117.31755, 40.07052, -276.33104, -65.08144, 59.46788, 60.81853,  
57                 19.97912, 47.13061, -71.74223, 5.18766, -76.9569, -9.07506, 3.9718 ),
58               c(53.18839, -97.95571, 249.06888, 52.84158, -61.07785, -30.62449,  
59                 0.18276, 37.84789, -220.10131, -20.13843, 7.1971, -25.77725,  105.98583),
60               c(-554.56612, -440.73448, -823.2569, -86.44144, 120.04677, -317.86967,  
61                 -140.17828, -116.4752, -556.48678, 4.49604, -22.66227, -66.98681,  -81.5922),
62               c(-15.57648, -261.64554, -168.64094, 307.14066, -150.18689, -289.77369,  
63                 44.92952, -47.7352, -14.03064, -111.88884, -0.19232, -70.95518,  91.73016))
64
65 trials <- 10
66 diff <- matrix(NA, trials, 13)
67
68 for (seed in 1:trials) {
69   set.seed(seed)
70   
71   m.mat <- mxMatrix(name="mean", nrow=1, ncol=6, values=runif(6, -1, 1), free=TRUE)
72   var <- runif(2, .1, 3)
73   cov1 <- runif(1, -min(abs(var)), min(abs(var)))
74   cov <- diag(runif(6, .1, 3))
75   cov[1:2,1:2] <- matrix(c(var[1],cov1,cov1,var[2]), nrow=2)
76   
77   cov.mat <- mxMatrix(name="cov", nrow=6, ncol=6, values=cov)
78   cov.mat@free <- cov.mat@values != 0
79   cov.mat@labels[1:2,1:2] <- matrix(c("v1","c12","c12","v2"), nrow=2)
80
81   m1 <- mxModel(model="latent",
82                 ip.mat, m.mat, cov.mat,
83                 mxData(observed=data, type="raw"),
84                 mxExpectationBA81(mean="mean", cov="cov",
85                                   ItemSpec=items, design=design,
86                                   ItemParam="ItemParam", qwidth=5, qpoints=21),
87                 mxFitFunctionML(),
88                 mxComputeSequence(steps=list(
89                   mxComputeOnce('expectation'),
90                   mxComputeOnce('fitfunction', gradient=TRUE))))
91   
92   objective1 <- function(x) {
93     m.mat@values[1:6] <- x[1:6]
94     diag(cov.mat@values) <- x[7:12]
95     cov.mat@values[1,2] <- x[13]
96     cov.mat@values[2,1] <- x[13]
97     m1.probe <- mxModel(m1, m.mat, cov.mat,
98                         mxComputeSequence(steps=list(
99                           mxComputeOnce('expectation'),
100                           mxComputeOnce('fitfunction', fit=TRUE))))
101     m1.probe <- mxRun(m1.probe, silent=TRUE)
102     got <- m1.probe@output$minimum
103     #  print(got)
104     got
105   }
106   
107   probe.pt <- c(m.mat@values, diag(cov.mat@values), cov.mat@values[2,1])
108 #  print(probe.pt - true.pt)
109   
110   deriv <- c()
111   if (0) {
112     require(numDeriv)
113     deriv <- grad(objective1, probe.pt, method="simple")
114     cat(deparse(round(deriv,5)), fill=TRUE)
115   } else {
116     deriv <- cache[[seed]]
117   }
118
119   m1 <- mxRun(m1, silent=TRUE)
120   #  print(probe.pt)
121   aGrad <- c(m1@output$gradient[-8], m1@output$gradient[8])
122   
123   diff[seed,] <- aGrad - deriv
124 }
125 #apply(abs(diff), 2, median)
126 omxCheckCloseEnough(apply(abs(diff), 2, median), rep(0,13), 2.6)