-
Notifications
You must be signed in to change notification settings - Fork 2
/
ora.R
40 lines (30 loc) · 941 Bytes
/
ora.R
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
### Prepare data
x <- unique(Orange$age)
y <- xtabs(circumference~as.character(Tree)+age, Orange)
y <- as.matrix(as.data.frame(as.matrix(unclass(y))))
dimnames(y) <- NULL
data <- list(x=x, y=y)
parameters <- list(phi1=200, phi2=800, phi3=400, b=rep(0,nrow(y)),
logSigma=0, logSigmaB=0)
### Compile model
require(TMB)
compile("ora.cpp")
dyn.load(dynlib("ora"))
### Run model
model <- MakeADFun(data, parameters, DLL="ora", random="b")
time <- system.time(fit <- nlminb(model$par, model$fn, model$gr))
best <- model$env$last.par.best
rep <- sdreport(model)
print(best)
print(rep)
print(summary(rep))
### Plot
phi1 <- best[["phi1"]]
phi2 <- best[["phi2"]]
phi3 <- best[["phi3"]]
b <- best[names(best)=="b"]
yfit <- matrix(NA_real_, nrow=5, ncol=7)
for(i in 1:5)
yfit[i,] <- (phi1 + b[i]) / (1+exp(-(x-phi2)/phi3))
matplot(x, t(y), xlim=c(0,1.1*max(x)), ylim=c(0,1.1*max(y)), ylab="y")
matlines(x, t(yfit), lty=1)