Demonstrate how to compute the Fisher information matrix
[openmx:openmx.git] / models / nightly / ifa-prob-gradient.R
1 # Demonstrate how to convert the derivative from log likelihood
2 # units to probability units.
3
4 library(rpf)
5 library(OpenMx)
6 library(numDeriv)
7
8 numItems <- 5
9 spec <- list()
10 spec[1:numItems] <- rpf.grm()
11 param <- sapply(spec, rpf.rparam)
12
13 pat <- data.frame(t(round(runif(numItems))))
14 for (c in colnames(pat)) {
15   pat[[c]] <- ordered(pat[[c]], levels=0:1)
16 }
17
18 if (0) {
19   data.fm <- sapply(rpf.sample(1000, spec, param), unclass)-1
20   write.table(data.fm, file="info.csv", quote=FALSE, row.names=FALSE, col.names=FALSE)
21 }
22
23 ip.mat <- mxMatrix(name="ItemParam", nrow=dim(param)[1], ncol=numItems,
24                    values=param, free=FALSE)
25
26 m.mat <- mxMatrix(name="mean", nrow=1, ncol=1, values=0, free=FALSE)
27 cov.mat <- mxMatrix(name="cov", nrow=1, ncol=1, values=1, free=FALSE)
28
29 objective <- function(x) {
30   ip.mat@values[,] <- x
31   m1 <- mxModel(model="latent",
32                 ip.mat, m.mat, cov.mat,
33                 mxData(observed=pat, type="raw"),
34                 mxExpectationBA81(mean="mean", cov="cov",
35                                   ItemSpec=spec,
36                                   ItemParam="ItemParam"),
37                 mxFitFunctionML(),
38                 mxComputeSequence(steps=list(
39                   mxComputeOnce('expectation'),
40                   mxComputeOnce('fitfunction', fit=TRUE))))
41   m1.fit <- mxRun(m1, silent=TRUE)
42   lik <- exp(m1.fit@output$minimum/-2)
43   lik
44 }
45
46 deriv <- grad(objective, c(param), method="simple")
47
48 m1 <- mxModel(model="latent",
49               ip.mat, m.mat, cov.mat,
50               mxData(observed=pat, type="raw"),
51               mxExpectationBA81(mean="mean", cov="cov",
52                                 ItemSpec=spec,
53                                 ItemParam="ItemParam"),
54               mxFitFunctionML(),
55               mxComputeSequence(steps=list(
56                 mxComputeOnce('expectation', context="EM"),
57                 mxComputeOnce('fitfunction', gradient=TRUE))))
58 m1@matrices$ItemParam@free[,] <- TRUE
59 m1.fit <- mxRun(m1, silent=TRUE)
60
61 ref <- m1.fit@expectation@patternLikelihood
62 grad <- exp(ref) - exp(ref + m1.fit@output$gradient)
63 omxCheckCloseEnough(grad, deriv, .01)   # poor accuracy
64 if (0) {
65   max(abs(grad- deriv))
66 }