ifa: snapshot with irt-2d working
[openmx:openmx.git] / models / failing / irt-2d.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
9 set.seed(7)
10
11 numItems <- 5
12 numPeople <- 500
13 maxDim <- 2
14
15 items <- vector("list", numItems)
16 correct <- vector("list", numItems)
17 for (ix in 1:numItems) {
18         items[[ix]] <- rpf.drm(dimensions=min(ix,maxDim),
19                                multidimensional=TRUE)
20         correct[[ix]] <- rpf.rparam(items[[ix]])
21         correct[[ix]][[4]] <- 0   # no guessing, for now
22 }
23 correct[[1]][[3]] <- 0   # no guessing, for now
24 correct.mat <- simplify2array(correct)
25
26 maxParam <- max(vapply(items, function(i) i@numParam, 0))
27 maxOutcomes <- max(vapply(items, function(i) i@numOutcomes, 0))
28
29 design <- matrix(c(1, 1,1,1,2,
30                   NA,2,2,2,1), byrow=TRUE, nrow=2)
31
32 ability <- array(rnorm(numPeople * 2), dim=c(numPeople, 2))
33 data <- rpf.sample(ability, items, correct, design)
34
35 spec <- mxMatrix(name="ItemSpec", nrow=6, ncol=numItems,
36          values=sapply(items, function(m) slot(m,'spec')),
37          free=FALSE, byrow=TRUE)
38
39 design <- mxMatrix(name="Design", nrow=maxDim, ncol=numItems,
40                   values=design)
41
42 ip.mat <- mxMatrix(name="ItemParam", nrow=maxParam, ncol=numItems,
43                    values=c(1, 1.4, 0, 0),
44                   lbound=c(1e-6, -1e6, 0, 0,
45                     rep(c(1e-6, 1e-6, -1e6, 0), 4)),
46          free=c(TRUE, TRUE, FALSE, FALSE,
47           rep(c(TRUE, TRUE, TRUE, FALSE), 4)))
48
49 m1 <- mxModel(model="2dim",
50           spec, design,
51           ip.mat,
52           mxData(observed=data, type="raw"),
53           mxExpectationBA81(
54              ItemSpec="ItemSpec",
55              Design="Design",
56              ItemParam="ItemParam",
57             GHpoints=21),
58           mxFitFunctionBA81())
59
60 m1 <- mxOption(m1, "Analytic Gradients", 'no')
61 if (1) {
62         m1 <- mxOption(m1, "Analytic Gradients", 'yes')
63         m1 <- mxOption(m1, "Verify level", '-1')
64 }
65 m1 <- mxOption(m1, "Function precision", '1.0E-5')
66 m1 <- mxOption(m1, "Calculate Hessian", "No")
67 m1 <- mxOption(m1, "Standard Errors", "No")
68
69 if (1) {
70
71         m1 <- mxRun(m1, silent=TRUE)
72
73         #print(m1@matrices$ItemParam@values)
74         #print(correct.mat)
75
76         omxCheckCloseEnough(cor(c(m1@matrices$ItemParam@values),
77                                 c(correct.mat)), .941, .01)
78         max.se <- max(m1@output$ability[,c(2,4)])
79         omxCheckCloseEnough(m1@output$ability[,c(1,3)], ability, max.se*2.5)
80         omxCheckCloseEnough(.576, cor(c(m1@output$ability[,c(1,3)]), c(ability)), .01)
81
82 } else { ################
83
84         myrpf <- function(iparam, spec, where) {
85                 ret <- array(dim=c(maxOutcomes, numItems))
86                 for (ix in 1:numItems) {
87                         i <- items[[ix]]
88                         param <- iparam[1:i@numParam,ix]
89                         out <- rpf.logprob(i, param, t(where[1:i@dimensions,ix]))
90                         ret[1:i@numOutcomes,ix] <- out
91                 }
92                 ret
93         }
94
95         m2 <- mxModel(name="R", m1,
96                       mxExpectationBA81(
97                         ItemSpec="ItemSpec",
98                         Design="Design",
99                         ItemParam="ItemParam",
100                         ItemPrior="prior",
101                         RPF=myrpf))
102
103         m2 <- mxRun(m2, silent=TRUE) # 23min, ugh
104         omxCheckCloseEnough(cor(c(m2@matrices$ItemParam@values),
105                                 c(correct.mat)), .934, .01)
106 }
107
108 q()