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