Fill in default upper & lower bounds
[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 spec <- mxMatrix(name="ItemSpec", nrow=3, ncol=numItems,
24          values=sapply(items, function(m) slot(m,'spec')),
25          free=FALSE, byrow=TRUE)
26
27 ip.mat <- mxMatrix(name="itemParam", nrow=4, ncol=numItems,
28                    values=c(1,0,0, 1),
29                    free=c(TRUE, TRUE, FALSE, FALSE))
30
31 eip.mat <- mxAlgebra(itemParam, name="EItemParam")
32
33 m.mat <- mxMatrix(name="mean", nrow=1, ncol=1, values=0, free=FALSE)
34 cov.mat <- mxMatrix(name="cov", nrow=1, ncol=1, values=1, free=FALSE)
35
36 m2 <- mxModel(model="drm1", ip.mat, spec, m.mat, cov.mat, eip.mat,
37               mxData(observed=data, type="raw"),
38               mxExpectationBA81(
39                 ItemSpec="ItemSpec", EItemParam="EItemParam",
40                 mean="mean", cov="cov",
41                 qpoints=31,
42                 scores="full"),
43               mxFitFunctionBA81(ItemParam="itemParam"),
44               mxComputeIterate(steps=list(
45                                  mxComputeOnce("EItemParam"),
46                                  mxComputeOnce('expectation', context='E'),
47                                    mxComputeNewtonRaphson(free.set='itemParam'),
48 #                                mxComputeGradientDescent(free.set='param'),
49                                  mxComputeOnce('expectation', context='M'),
50                                  mxComputeOnce('fitfunction')
51                                  )))
52
53         m2 <- mxOption(m2, "Analytic Gradients", 'Yes')
54         m2 <- mxOption(m2, "Verify level", '-1')
55 m2 <- mxOption(m2, "Function precision", '1.0E-5')
56 m2 <- mxRun(m2)
57
58 #print(m2@matrices$itemParam@values)
59 #print(correct.mat)
60 got <- cor(c(m2@matrices$itemParam@values[1:2,]),
61            c(correct.mat[1:2,]))
62 omxCheckCloseEnough(got, .988, .01)
63 ability <- scale(ability)
64 scores <- m2@expectation@scores.out
65 omxCheckCloseEnough(scores[,1], as.vector(ability), 3.5*max(scores[,2]))
66 omxCheckCloseEnough(cor(c(scores[,1]), ability), .737, .01)