Gradients for the latent distribution
[openmx:openmx.git] / models / passing / ifa-latent.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(15)
11 correct.deriv <- c(484.533329156793, 892.742069416003, 156.534557121965, 283.449559864665,  9.29526394590284)
12 correct.bifactor <- c(-48.70145, 280.13851, 23.39651, 109.05235, 31.38746, 209.36438,
13                       172.0925, 124.71525, 164.80765, 131.46529)
14
15 numItems <- 20
16 numPeople <- 1000
17
18 items <- list()
19 correct <- list()
20 for (ix in 1:numItems) {
21         items[[ix]] <- rpf.grm(factors=2)
22         correct[[ix]] <- rpf.rparam(items[[ix]])
23 }
24 correct.mat <- simplify2array(correct)
25
26 maxParam <- max(vapply(items, function(i) rpf.numParam(i), 0))
27
28 true.mean <- c(.2,-.1)
29 true.cov <- matrix(c(.5, -.2, -.2, .5), nrow=2)
30 ability <- rmvnorm(numPeople, mean=true.mean, sigma=true.cov)
31 data <- rpf.sample(t(ability), items, correct.mat)
32
33 ip.mat <- mxMatrix(name="ItemParam", nrow=maxParam, ncol=numItems,
34                    values=correct.mat, free=FALSE)
35
36 eip.mat <- mxAlgebra(ItemParam, name="EItemParam")
37
38 m.mat <- mxMatrix(name="mean", nrow=1, ncol=2, values=c(.5, .5), free=TRUE)
39 cov.mat <- mxMatrix(name="cov", nrow=2, ncol=2, values=matrix(c(1,.2,.2,1), nrow=2),
40                     free=TRUE, labels=c("v1","c12","c12","v2"))
41
42 m1 <- mxModel(model="latent",
43               ip.mat, m.mat, cov.mat, eip.mat,
44               mxData(observed=data, type="raw"),
45               mxExpectationBA81(mean="mean", cov="cov",
46                                 ItemSpec=items,
47                                 EItemParam="EItemParam"),
48               mxFitFunctionBA81(ItemParam="ItemParam"),
49               mxComputeSequence(steps=list(
50                 mxComputeOnce("EItemParam"),
51                 mxComputeOnce('expectation'),
52                 mxComputeOnce('fitfunction', gradient=TRUE))))
53
54 if (1) {
55 #  m1 <- mxOption(m1, "Number of Threads", 1)
56 #  m1 <- mxOption(m1, "No Sort Data", 'latent')
57   m1 <- mxRun(m1, silent=TRUE)
58   omxCheckCloseEnough(m1@output$gradient, correct.deriv, 1e-3)
59 }
60
61 if (0) {
62   m1 <- mxModel(m1,
63                 mxComputeSequence(steps=list(
64                   mxComputeOnce("EItemParam"),
65                   mxComputeOnce('expectation'),
66                   mxComputeGradientDescent(useGradient=TRUE))))
67 #  m1 <- mxOption(m1, "Number of Threads", 1)
68   m1 <- mxOption(m1, "Verify level", '-1')
69   m1 <- mxOption(m1, "Function precision", '1.0E-5')
70   m1 <- mxRun(m1, silent=TRUE)
71   omxCheckCloseEnough(c(m1@matrices$mean@values), true.mean, .02)
72   omxCheckCloseEnough(m1@matrices$cov@values, true.cov, .1)
73 }
74
75 objective1 <- function(x) {
76         m.mat@values[1:2] <- x[1:2]
77         cov.mat@values[1,1] <- x[3]
78         cov.mat@values[1,2] <- x[4]
79         cov.mat@values[2,1] <- x[4]
80         cov.mat@values[2,2] <- x[5]
81         m1 <- mxModel(m1, m.mat, cov.mat,
82                       mxComputeSequence(steps=list(
83                                           mxComputeOnce("EItemParam"),
84                                           mxComputeOnce('expectation'),
85                                           mxComputeOnce('fitfunction'))))
86         m1 <- mxRun(m1, silent=TRUE)
87         got <- m1@output$minimum
88 #  print(got)
89   got
90 }
91
92 if (0) {
93   require(numDeriv)
94   deriv <- grad(objective1, c(m.mat@values, cov.mat@values[1:2,1], cov.mat@values[2,2]))
95   
96   omxCheckCloseEnough(c(deriv), correct.deriv, 1e-3)
97 }
98
99 if (1) {
100   m.mat <- mxMatrix(name="mean", nrow=1, ncol=5, values=rep(.1,5), free=TRUE)
101   cov.mat <- mxMatrix(name="cov", nrow=5, ncol=5, values=diag(seq(.81,1.41,length.out=5)),
102                       free=diag(5)==1)
103   
104   design <- matrix(c(rep(1,numItems),
105                      kronecker(2:5,rep(1,5))), byrow=TRUE, ncol=numItems)
106   m1 <- mxModel(m1, m.mat, cov.mat,
107                 mxMatrix(name="Design", nrow=dim(design)[1], ncol=numItems, values=design),
108                 mxExpectationBA81(mean="mean", cov="cov",
109                                   ItemSpec=items,
110                                   Design="Design",
111                                   EItemParam="EItemParam"))
112 #  m1 <- mxOption(m1, "Number of Threads", 1)
113   m1 <- mxRun(m1, silent=TRUE)
114   omxCheckCloseEnough(m1@output$gradient[1:5], correct.bifactor[1:5], 1e-3)  #means
115   omxCheckCloseEnough(m1@output$gradient[6:10], correct.bifactor[6:10], 1e-2) #vars
116 }
117
118 objective2 <- function(x) {
119   m.mat@values[1:5] <- x[1:5]
120   cov.mat@values <- diag(x[6:10])
121   m1 <- mxModel(m1, m.mat, cov.mat,
122                 mxComputeSequence(steps=list(
123                   mxComputeOnce("EItemParam"),
124                   mxComputeOnce('expectation'),
125                   mxComputeOnce('fitfunction'))))
126   m1 <- mxRun(m1, silent=TRUE)
127   got <- m1@output$minimum
128   #  print(got)
129   got
130 }
131
132 if (0) {
133   require(numDeriv)
134   deriv <- grad(objective2, c(m.mat@values, diag(cov.mat@values)))
135   cat(deparse(round(deriv, 5)))
136   omxCheckCloseEnough(c(deriv), correct.deriv, 1e-3)
137 }