Check expectation@patternLikelihood in log units
[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",ItemParam="itemParam",
36                 mean="mean", cov="cov", qpoints=31),
37               mxFitFunctionML(),
38               mxComputeOnce('expectation', context='EM'))
39 m2 <- mxRun(m2)
40 omxCheckCloseEnough(sum(m2@expectation@patternLikelihood), 0.3846289, .01)
41 omxCheckCloseEnough(fivenum(log(m2@expectation@patternLikelihood)),
42                     c(-7.5454472, -7.3950031, -7.3950031, -6.9391761, -3.5411989), .001)
43 omxCheckCloseEnough(sum(m2@expectation@em.expected), 5000, .01)
44 omxCheckCloseEnough(fivenum(m2@expectation@em.expected),
45                     c(0, 5.86e-05, 0.0687802, 7.1582354, 74.1583248), .01)
46
47 m2 <- mxModel(m2,
48               mxComputeIterate(steps=list(
49                                  mxComputeOnce('expectation', context='EM'),
50                                  mxComputeOnce('fitfunction', gradient=TRUE, hessian=TRUE)
51                                  )))
52 m2 <- mxRun(m2)
53 omxCheckCloseEnough(m2@fitfunction@result, 3221.826, .01)
54 omxCheckCloseEnough(fivenum(m2@output$gradient), c(-128.034, -8.294, 10.7, 25.814, 107.966), .01)
55 omxCheckCloseEnough(fivenum(m2@output$hessian[m2@output$hessian != 0]),
56                     c(6.559, 6.559, 32.976, 83.554, 107.714), .01)
57
58 m2 <- mxModel(m2,
59               mxData(observed=data, type="raw"),  # got sorted, add it again unsorted
60               mxExpectationBA81(
61                 ItemSpec=items, EItemParam="EItemParam",ItemParam="itemParam",
62                 mean="mean", cov="cov",
63                 qpoints=31,
64                 scores="full"),
65               mxComputeIterate(steps=list(
66                                  mxComputeOnce('expectation', context='EM'),
67                                  mxComputeNewtonRaphson(free.set='itemParam'),
68                                  mxComputeOnce('expectation'),
69                                  mxComputeOnce('fitfunction', free.set=c("mean","cov"))
70                                  )))
71
72         m2 <- mxOption(m2, "Analytic Gradients", 'Yes')
73         m2 <- mxOption(m2, "Verify level", '-1')
74 m2 <- mxOption(m2, "Function precision", '1.0E-5')
75 m2 <- mxRun(m2)
76
77 #print(m2@matrices$itemParam@values)
78 #print(correct.mat)
79 omxCheckCloseEnough(m2@fitfunction@result, 6216.272, .01)
80 got <- cor(c(m2@matrices$itemParam@values[1:2,]),
81            c(correct.mat[1:2,]))
82 omxCheckCloseEnough(got, .988, .01)
83 scores <- m2@expectation@scores.out
84 omxCheckCloseEnough(scores[1:5,1], c(0.6783773, 0.2848123, -0.3438632, -0.1026575, -1.0820213), .001)
85 omxCheckCloseEnough(scores[1:5,2], c(0.6769653, 0.6667262, 0.6629124, 0.6624804, 0.6796952), 1e-5)
86 omxCheckCloseEnough(scores[,1], as.vector(ability), 3.5*max(scores[,2]))
87 omxCheckCloseEnough(cor(c(scores[,1]), ability), .737, .01)