Update for evolving OpenMx API
[rpf:rpf.git] / util / test-dLL.R
1 library(testthat)
2 library(OpenMx)  # an unrelease version of OpenMx is required for this test
3 library(rpf)
4 library(numDeriv)
5 #options(error = utils::recover)
6
7 context("dLL")
8
9 unpackHession <- function(deriv, np) {
10   hess <- matrix(NA, nrow=np, ncol=np)
11   dx <- np+1
12   for (hr in 1:np) {
13     hess[1:hr,hr] <- hess[hr,1:hr] <- deriv$D[dx:(dx+hr-1)]
14     dx <- dx + hr
15   }
16   hess
17 }
18
19 myseed <- as.integer(runif(1) * 1e7)
20 #print(paste("set.seed =",myseed))
21 set.seed(myseed)
22 #set.seed(3090526)
23
24 numItems <- 3
25 items <- list()
26 items[[1]] <- rpf.drm(factors=2)
27 items[[2]] <- rpf.grm(outcomes=3, factors=2)
28 T.a <- matrix(rnorm(9),3,3) + diag(3)
29 T.c <- matrix(rnorm(9),3,3) + diag(3)
30
31 items[[3]] <- rpf.nrm(outcomes=4, factors=2,
32                       T.a=T.a, T.c=T.c)
33
34 # If not all outcomes are represented then lots of warnings result.
35 # This is not a cause for concern.
36 data <- rpf.sample(200, items)
37
38 starting <- list(c(1.4, 1, 0, .1, .9),
39                  c(1.4, 1, -.5, -1),
40                  c(1.4,  1,  rep(0,6)))
41 starting.len <- max(vapply(starting, length, 0))
42
43 ip.mat <- mxMatrix(name="itemParam", nrow=starting.len, ncol=numItems,
44                    values=0, free=FALSE)
45
46 for (sx in 1:length(starting)) {
47   v <- starting[[sx]]
48   ip.mat@values[1:length(v),sx] <- v
49   ip.mat@free[1:length(v),sx] <- TRUE
50 }
51
52 Eip <- mxMatrix(name="EitemParam", nrow=dim(ip.mat@values)[1], ncol=numItems,
53                 values=ip.mat@values)
54
55 m.mat <- mxMatrix(name="mean", nrow=1, ncol=2, values=0, free=FALSE)
56 cov.mat <- mxMatrix(name="cov", nrow=2, ncol=2, values=diag(2), free=FALSE)
57 m2 <- mxModel(model="drm1", ip.mat, Eip, m.mat, cov.mat,
58               mxData(observed=data, type="raw"),
59               mxExpectationBA81(mean="mean", cov="cov",
60                                 ItemSpec=items,
61                                 ItemParam="itemParam", EItemParam="EitemParam",
62                 qpoints=13),
63               mxFitFunctionML(),
64               mxComputeSequence(steps=list(
65                                   mxComputeOnce('expectation', context='EM'),
66                                   mxComputeOnce('fitfunction', gradient=TRUE, hessian=TRUE)
67                                   )))
68
69 spoint <- list(c(1.4, 1, 0, .1, .9),
70                c(1.4, 1, .5, -.5),
71                c(0.89,  0.33,  0.05,  0.18,  0.07, -2.03,  0.17,  0.26))
72 spoint.len <- vapply(spoint, length, 0)
73
74 for (ii in 1:numItems) {
75   np <- length(spoint[[ii]])
76   m2@matrices$itemParam@values <- Eip@values
77   
78   m2@matrices$itemParam@values[1:np,ii] <- spoint[[ii]]
79   m2 <- mxRun(m2, silent=TRUE)
80   
81   offset <- 1
82   if (ii > 1) offset <- sum(c(1,spoint.len[1:(ii-1)]))
83   offset.i <- seq(offset,offset+np-1)
84
85   grad1 <- m2@output$gradient[offset.i]
86   names(grad1) <- NULL
87   hess <- m2@output$hessian[offset.i, offset.i]
88   for (cx in 2:dim(hess)[1]) for (rx in 1:cx) hess[rx,cx] <- hess[cx,rx]
89   
90   deriv <- genD(function(param) {
91     np <- length(param)
92     m2@matrices$itemParam@values[1:np,ii] <- param
93     m2 <- mxRun(m2, silent=TRUE)
94     m2@output$minimum
95   }, spoint[[ii]], method.args=list(eps=0.01, d=0.01, r=2))
96
97   emp.hess <- unpackHession(deriv, np)
98   emp.hess[is.na(emp.hess)] <- 0
99
100   if (0) {
101     print(paste("Item", ii))
102     print(grad1)
103     print(deriv$D[1:np])
104     cat("T.a=",deparse(T.a),"\n")
105     cat("T.c=",deparse(T.c),"\n")
106     cat("an=",deparse(hess),"\n")
107     cat("emp=",deparse(emp.hess),"\n")
108     print(round(hess - emp.hess, 2))
109   }
110   
111   expect_equal(deriv$D[1:np], grad1, tolerance=1e-6)
112   expect_equal(emp.hess, hess, tolerance=1e-4)
113 }
114 #warnings()