1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
|
## ---- include = FALSE---------------------------------------------------------
is_CRAN <- !identical(Sys.getenv("NOT_CRAN"), "true")
if (!is_CRAN) {
options(mc.cores = parallel::detectCores())
} else {
knitr::opts_chunk$set(eval = FALSE)
knitr::knit_hooks$set(evaluate.inline = function(x, envir) x)
}
## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
## -----------------------------------------------------------------------------
library(OpenMx)
aPlan <- mxComputeSequence(list( #analytic
mxComputeOnce('fitfunction', c('gradient')),
mxComputeReportDeriv()))
nPlan <- mxComputeSequence(list( #numerical
mxComputeNumericDeriv(analytic = FALSE, hessian=FALSE, checkGradient = FALSE),
mxComputeReportDeriv()))
## -----------------------------------------------------------------------------
mat1 <- mxMatrix("Full", rnorm(1), free=TRUE, nrow=1, ncol=1, labels="m1", name="mat1")
obj <- mxAlgebra(-.5 * (log(2*pi) + log(2) + (mat1[1,1])^2/2), name = "obj")
grad <- mxAlgebra(-(mat1[1,1]) + 2, name = "grad", dimnames=list("m1", NULL))
mv1 <- mxModel("mv1", mat1, obj, grad,
mxFitFunctionAlgebra("obj", gradient="grad"))
## -----------------------------------------------------------------------------
nu <- mxRun(mxModel(mv1, nPlan), silent = TRUE)
an <- mxRun(mxModel(mv1, aPlan), silent = TRUE)
cbind(numerical=nu$output$gradient, analytic=an$output$gradient)
## -----------------------------------------------------------------------------
p1 <- runif(2, -10,10)
mv1 <- omxSetParameters(mv1, labels = 'm1', values=p1)
nu <- mxRun(mxModel(mv1, nPlan), silent = TRUE)
an <- mxRun(mxModel(mv1, aPlan), silent = TRUE)
cbind(numerical=nu$output$gradient, analytic=an$output$gradient)
## -----------------------------------------------------------------------------
grad <- mxAlgebra(-(mat1[1,1])/2, name = "grad", dimnames=list("m1", NULL))
mv2 <- mxModel(mv1, grad)
## -----------------------------------------------------------------------------
nu <- mxRun(mxModel(mv2, nPlan), silent = TRUE)
an <- mxRun(mxModel(mv2, aPlan), silent = TRUE)
cbind(numerical=nu$output$gradient, analytic=an$output$gradient)
## -----------------------------------------------------------------------------
mv2 <- omxSetParameters(mv2, labels = 'm1', values=rnorm(1))
nu <- mxRun(mxModel(mv2, nPlan), silent = TRUE)
an <- mxRun(mxModel(mv2, aPlan), silent = TRUE)
cbind(numerical=nu$output$gradient, analytic=an$output$gradient)
|
