Gradients for the latent distribution
[openmx:openmx.git] / models / passing / ifa-drm1.R
1 options(error = utils::recover)
2 library(OpenMx)
3 library(rpf)
4
5 set.seed(9)
6
7 numItems <- 10
8 i1 <- rpf.drm(multidimensional=TRUE)
9 items <- vector("list", numItems)
10 correct <- vector("list", numItems)
11 for (ix in 1:numItems) {
12   items[[ix]] <- i1
13   correct[[ix]] <- rpf.rparam(i1)
14   correct[[ix]][3] <- 0
15   correct[[ix]][4] <- 1
16 }
17 correct.mat <- simplify2array(correct)
18 correct.mat[2,] <- correct.mat[2,] * -correct.mat[1,]
19
20 ability <- rnorm(500)
21 data <- rpf.sample(ability, items, correct.mat)
22
23 ip.mat <- mxMatrix(name="itemParam", nrow=4, ncol=numItems,
24                    values=c(1,0,0, 1),
25                    free=c(TRUE, TRUE, FALSE, FALSE))
26
27 eip.mat <- mxAlgebra(itemParam, name="EItemParam")
28
29 m.mat <- mxMatrix(name="mean", nrow=1, ncol=1, values=0, free=FALSE)
30 cov.mat <- mxMatrix(name="cov", nrow=1, ncol=1, values=1, free=FALSE)
31
32 m2 <- mxModel(model="drm1", ip.mat, m.mat, cov.mat, eip.mat,
33               mxData(observed=data, type="raw"),
34               mxExpectationBA81(
35                 ItemSpec=items, EItemParam="EItemParam",
36                 mean="mean", cov="cov",
37                 qpoints=31,
38                 scores="full"),
39               mxFitFunctionBA81(ItemParam="itemParam"),
40               mxComputeIterate(steps=list(
41                                  mxComputeOnce("EItemParam"),
42                                  mxComputeOnce('expectation', context='EM'),
43                                  mxComputeNewtonRaphson(free.set='itemParam'),
44                                  mxComputeOnce('expectation'),
45                                  mxComputeOnce('fitfunction')
46                                  )))
47
48         m2 <- mxOption(m2, "Analytic Gradients", 'Yes')
49         m2 <- mxOption(m2, "Verify level", '-1')
50 m2 <- mxOption(m2, "Function precision", '1.0E-5')
51 m2 <- mxRun(m2)
52
53 #print(m2@matrices$itemParam@values)
54 #print(correct.mat)
55 got <- cor(c(m2@matrices$itemParam@values[1:2,]),
56            c(correct.mat[1:2,]))
57 omxCheckCloseEnough(got, .988, .01)
58 ability <- scale(ability)
59 scores <- m2@expectation@scores.out
60 omxCheckCloseEnough(scores[,1], as.vector(ability), 3.5*max(scores[,2]))
61 omxCheckCloseEnough(cor(c(scores[,1]), ability), .737, .01)