Compute condition number of information matrix
[openmx:openmx.git] / models / nightly / ifa-dLL.R
1 #options(digits=20)
2 library(OpenMx)
3 library(rpf)
4 library(numDeriv)
5 #options(error = utils::recover)
6
7 unpackHession <- function(deriv, np) {
8   hess <- matrix(NA, nrow=np, ncol=np)
9   dx <- np+1
10   for (hr in 1:np) {
11     hess[1:hr,hr] <- hess[hr,1:hr] <- deriv[dx:(dx+hr-1)]
12     dx <- dx + hr
13   }
14   hess
15 }
16
17 #myseed <- as.integer(runif(1) * 1e7)
18 #print(paste("set.seed =",myseed))
19 #set.seed(myseed)
20 set.seed(1)
21
22 items <- list()
23 items[[1]] <- rpf.drm(factors=2)
24 items[[2]] <- rpf.grm(outcomes=5, factors=2)
25 items[[3]] <- rpf.nrm(outcomes=4, factors=2,
26                       T.a=diag(3), T.c=diag(3))
27 numItems <- length(items)
28
29 params <- lapply(items, rpf.rparam)
30 data <- rpf.sample(5000, items, params)
31
32 starting <- list(c(1.4, 1, 0, .1, .9),
33                  c(1.4, 1, seq(2,-2, length.out=4)),
34                  c(1.4,  1,  rep(0,6)))
35 starting.len <- max(vapply(starting, length, 0))
36
37 ip.mat <- mxMatrix(name="itemParam", nrow=starting.len, ncol=numItems,
38                    values=0, free=FALSE)
39
40 for (sx in 1:length(starting)) {
41   v <- starting[[sx]]
42   ip.mat@values[1:length(v),sx] <- v
43   ip.mat@free[1:length(v),sx] <- TRUE
44 }
45 starting.free <- ip.mat@free
46 starting.values <- ip.mat@values
47
48 m.mat <- mxMatrix(name="mean", nrow=1, ncol=2, values=0, free=FALSE)
49 cov.mat <- mxMatrix(name="cov", nrow=2, ncol=2, values=diag(2), free=FALSE)
50 m2 <- mxModel(model="drm1", ip.mat, m.mat, cov.mat,
51               mxData(observed=data, type="raw"),
52               mxExpectationBA81(mean="mean", cov="cov",
53                                 ItemSpec=items,
54                                 ItemParam="itemParam",
55                                 EItemParam=starting.values,
56                 qwidth=5, qpoints=21),
57               mxFitFunctionML())
58
59 if (0) {   # enable to generate answer file
60   samples.per.item <- 100
61   ans <- list()
62   for (ii in 1:numItems) {  # 
63     spi <- 0
64     while (spi < samples.per.item) {
65       m2@matrices$itemParam@values <- starting.values
66       m2@matrices$itemParam@free[,] <- FALSE
67       
68       spoint <- rpf.rparam(items[[ii]])
69       # exclude GRM with close adjacent intercepts, too much curvature for numDeriv
70       if (ii==2 && min(abs(diff(spoint[3:6]/max(spoint[1:2])))) < .3) next
71
72       np <- length(spoint)
73       m2@matrices$itemParam@values[1:np,ii] <- spoint
74       
75       deriv <- genD(function(param) {
76         np <- length(param)
77         m2@matrices$itemParam@values[1:np,ii] <- param
78         lModel <- mxModel(m2,
79                           mxComputeSequence(steps=list(
80                             mxComputeOnce('expectation', context='EM'),
81                             mxComputeOnce('fitfunction', fit=TRUE)
82                           )))
83         fit <- mxRun(lModel, silent=TRUE)
84         fit@output$minimum
85       }, spoint, method.args=list(d=.01, r=2))
86       
87       # Our gradients are too flat for the default higher precision settings.
88       
89       if (any(is.na(deriv$D))) next
90
91       print(c(ii, spoint))
92       ans[[length(ans)+1]] <- c(ii, spoint, deriv$D)
93       spi <- spi + 1
94     }
95   }
96   
97   ans.len <- max(sapply(ans, length))
98   ans.padded <- lapply(ans, function (elem) elem <- c(elem, rep(NA, ans.len-length(elem))))
99   write.table(t(simplify2array(ans.padded)), file="data/dLL.csv", row.names=FALSE, col.names=FALSE)
100 }
101
102 ans <- suppressWarnings(try(read.table("models/nightly/data/dLL.csv"), silent=TRUE))
103 if (is(ans, "try-error")) ans <- read.table("data/dLL.csv")
104
105 m2 <- mxModel(m2,
106               mxComputeSequence(steps=list(
107                 mxComputeOnce('expectation', context='EM'),
108                 mxComputeOnce('fitfunction', gradient=TRUE, hessian=TRUE)
109               )))
110
111 if (1) {  # enable to examine the RMSE by item model
112   for (ix in 1:numItems) {
113     np <- rpf.numParam(items[[ix]])
114     diff.grad <- rep(0, np)
115     diff.hess <- matrix(0, np, np)
116     diff.count <- 0
117     skip <- c()
118     
119     for (tx in 1:dim(ans)[1]) {
120       ii <- ans[tx,1]
121       if (ii != ix) next
122       
123       m2@matrices$itemParam@values <- starting.values
124       m2@matrices$itemParam@free[,] <- FALSE
125       
126       spoint <- ans[tx,2:(np+1)]
127       m2@matrices$itemParam@values[1:np,ii] <- simplify2array(spoint)
128       m2@matrices$itemParam@free[,ii] <- starting.free[,ii]
129       
130       m2 <- mxRun(m2, silent=TRUE)
131       
132       grad1 <- m2@output$gradient
133       names(grad1) <- NULL
134       hess <- m2@output$hessian
135       
136       emp.grad <- simplify2array(ans[tx,(2+np):(1+2*np)])
137       emp.hess <- unpackHession(simplify2array(ans[tx, -1:-(1+np)]), np)
138       
139       if (any(abs(emp.hess - hess) > .01)) {
140         skip <- c(skip,tx)
141       }
142       diff.grad <- diff.grad + (emp.grad - grad1)^2
143       diff.hess <- diff.hess + (emp.hess - hess)^2
144       diff.count <- diff.count+1
145     }
146
147     if (diff.count > 0 && FALSE) {
148       if(length(skip)) print(skip)
149       diff.grad <- sqrt(diff.grad / diff.count)
150       diff.hess <- sqrt(diff.hess / diff.count)
151       print(diff.grad)
152       print(max(diff.grad))
153       print(diff.hess)
154       print(max(diff.hess))
155     }
156 #     print(max(diff.grad))
157 #     print(max(diff.hess))
158     # The poor accuracy here is probably due to numDeriv, not the
159     # math for analytic derivs.
160     omxCheckTrue(all(diff.grad < .003))
161     omxCheckTrue(all(diff.hess < .11))
162   }
163 }
164
165 if (0) {
166   kat <- c()
167   badest <-  c(6, 15, 40, 55, 67, 82, 84, 85, 87, 94)
168  # badest <- c(107, 114, 121, 126, 132, 138, 139, 164, 172, 177, 182, 199)
169       for (tx in badest) {
170     ii <- ans[tx,1]
171 #    print(params[[ii]])
172     np <- rpf.numParam(items[[ii]])
173     
174     m2@matrices$itemParam@values <- starting.values
175     m2@matrices$itemParam@free[,] <- FALSE
176     
177     spoint <- simplify2array(ans[tx,2:(np+1)])
178     print(spoint)
179       next;
180   
181     m2@matrices$itemParam@values[1:np,ii] <- spoint
182     m2@matrices$itemParam@free[,ii] <- starting.free[,ii]
183       
184     m2 <- mxRun(m2, silent=TRUE)
185       
186     grad1 <- m2@output$gradient
187     names(grad1) <- NULL
188     hess <- m2@output$hessian
189     
190     if (0) {
191       print(paste("Item", ii))
192       print(grad1)
193       print(deriv$D[1:np])
194       cat("T.a=",deparse(T.a),"\n")
195       cat("T.c=",deparse(T.c),"\n")
196       cat("an=",deparse(hess),"\n")
197       cat("emp=",deparse(emp.hess),"\n")
198       print(round(hess - emp.hess, 2))
199     }
200     
201     emp.grad <- simplify2array(ans[tx,(2+np):(1+2*np)])
202     emp.hess <- unpackHession(simplify2array(ans[tx, -1:-(1+np)]), np)
203
204 #    print(grad1)
205 #    print(grad1 - emp.grad)
206 #    print(emp.hess)
207 #    print(hess)
208 #    print(hess - emp.hess)
209     
210     evalLL <- function(param) {
211       np <- length(param)
212       m2@matrices$itemParam@values[1:np,ii] <- param
213       lModel <- mxModel(m2,
214                         mxComputeSequence(steps=list(
215                           mxComputeOnce('expectation', context='EM'),
216                           mxComputeOnce('fitfunction', fit=TRUE)
217                         )))
218       fit <- mxRun(lModel, silent=TRUE)
219       ll <- fit@output$minimum
220       ll
221     }
222     deriv <- genD(evalLL, spoint, method.args=list(d=.1, r=3))
223     print(round(hess - unpackHession(deriv$D, np), 2))
224     
225     if (0) {
226       grid <- expand.grid(x=seq(.07,-.05,-.01))
227       for (pp in grid$x) {
228         spoint[5] <- pp
229         grid$LL <- evalLL(spoint)
230       }
231     }
232   }
233 }