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