Move all fitfunction args to expectation
[openmx:openmx.git] / models / passing / ifa-2d.R
1 # echo 0 > /proc/self/coredump_filter  # normally 023
2 # R --vanilla --no-save -f models/failing/bock-aitkin-1981.R
3 # R -d gdb --vanilla --no-save -f models/failing/bock-aitkin-1981.R
4
5 #options(error = browser)
6 require(OpenMx)
7 require(rpf)
8 #library(mvtnorm)
9
10 set.seed(7)
11
12 numItems <- 5
13 numPeople <- 500
14 maxDim <- 2
15
16 items <- vector("list", numItems)
17 correct <- vector("list", numItems)
18 for (ix in 1:numItems) {
19         items[[ix]] <- rpf.drm(factors=min(ix,maxDim),
20                                multidimensional=TRUE)
21         correct[[ix]] <- rpf.rparam(items[[ix]])
22         if (ix>1) correct[[ix]][[4]] <- 0   # no guessing, for now
23 }
24 correct[[1]][[5]] <- 1   # make all vectors the same length
25 correct.mat <- simplify2array(correct)
26 correct.mat[5,] <- 1
27 correct.mat[4,1] <- 1
28 correct.mat[3,1] <- 0
29
30 maxParam <- max(vapply(items, function(i) rpf.numParam(i), 0))
31 maxOutcomes <- max(vapply(items, function(i) i@outcomes, 0))
32
33 design <- matrix(c(1, 1,1,1,2,
34                   NA,2,2,2,1), byrow=TRUE, nrow=2)
35
36 ability <- array(rnorm(numPeople * 2), dim=c(2, numPeople))
37 #cov <- matrix(c(1, .68, .68, 1), nrow=2)
38 #ability <- rmvnorm(numPeople, sigma=cov)
39 data <- rpf.sample(ability, items, correct.mat, design)
40 #write.csv(data, file="2d-new.csv")
41
42 ip.mat <- mxMatrix(name="ItemParam", nrow=maxParam, ncol=numItems,
43                    values=c(1, 1.4, 0, 0, 1),
44          free=c(TRUE, TRUE, FALSE, FALSE, FALSE,
45           rep(c(TRUE, TRUE, TRUE, FALSE, FALSE), 4)))
46 ip.mat@values[4,1] <- 1
47
48 eip.mat <- mxAlgebra(ItemParam, name="EItemParam")
49
50 m.mat <- mxMatrix(name="mean", nrow=1, ncol=2, values=0, free=FALSE)
51 cov.mat <- mxMatrix(name="cov", nrow=2, ncol=2, values=diag(2), free=FALSE)
52
53 m1 <- mxModel(model="2dim",
54           ip.mat, m.mat, cov.mat, eip.mat,
55           mxData(observed=data, type="raw"),
56           mxExpectationBA81(mean="mean", cov="cov",
57              ItemSpec=items,
58                             ItemParam="ItemParam", EItemParam="EItemParam",
59              design=design,
60             qpoints=29,
61             scores="full"),
62           mxFitFunctionBA81(),
63               mxComputeIterate(steps=list(
64                 mxComputeOnce("EItemParam"),
65                 mxComputeOnce('expectation', context='EM'),
66                                  mxComputeNewtonRaphson(free.set='ItemParam'),
67 #                                mxComputeGradientDescent(free.set='ItemParam'),
68                 mxComputeOnce('expectation'),
69                 mxComputeOnce('fitfunction'))))
70
71         m1 <- mxOption(m1, "Analytic Gradients", 'Yes')
72         m1 <- mxOption(m1, "Verify level", '-1')
73 m1 <- mxOption(m1, "Function precision", '1.0E-5')
74
75 m1 <- mxRun(m1, silent=TRUE)
76
77                                         #print(m1@matrices$ItemParam@values)
78                                         #print(correct.mat)
79 # sometimes found as low as .89, maybe solution is unstable
80 omxCheckCloseEnough(cor(c(m1@matrices$ItemParam@values),
81                         c(correct.mat)), .9, .06)
82 scores.out <- m1@expectation@scores.out
83 max.se <- max(scores.out[,3:4])
84 omxCheckCloseEnough(sum(abs(scores.out[,1:2] - t(ability)) < max.se) / (numPeople*2), .797, .02)
85 omxCheckCloseEnough(.672, cor(c(scores.out[,1:2]), c(t(ability))), .01)