-
Notifications
You must be signed in to change notification settings - Fork 1
/
simulate.R
81 lines (75 loc) · 1.62 KB
/
simulate.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
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
70
71
72
73
74
75
76
77
78
79
80
# simulate.R
beetle_sim <- function( state = c(100, 10, 10, 10, 10),
pars = c( 5,
0,
0.001,
0.0,
0.003,
3.8,
3.8+16.4,
3.8+16.4+5,
0.01,
0.004,
0.01,
3.8+8
),
dt = 14,
T = 200,
seed = 123 ){
out <- .C("_Z10beetle_simPiPdS0_S0_S_", as.integer(state), as.double(pars), as.double(dt), as.double(T), as.integer(seed) )
data <- read.table("beetle_sim.txt")
list(state = data, pars = out[[2]], dt = dt)
}
ensemble <- function( state = c(450, 0, 0, 0, 30),
initial = state,
pars = c( 5,
0,
0.001,
0.0,
0.003,
3.8,
3.8+16.4,
3.8+16.4+5,
0.01,
0.004,
0.01,
3.8+8
),
dt = 14,
seed = 123,
reps = 10
)
{
probs = double(5)
out <- .C( "_Z8ensemblePiS_PdS0_S_S_S0_S_",
as.integer(state),
as.integer(initial),
as.double(pars),
as.double(dt),
as.integer(seed),
as.integer(reps),
as.double(probs),
as.integer(5) )
list(probs=out[[7]], pars = out[[2]], dt = dt)
}
likelihood <- function(pars, X, dt, reps=100, cpus=2){
# compute likelihood
n <- length(X[,1])
starts <- X[1:n-1, ]
ends <- X[2:n, ]
seed <- 1e9*runif(n-1)
require(snowfall)
if (cpus > 1){
sfInit(parallel=TRUE, cpus=cpus)
} else {
sfInit()
}
sfLibrary(beetles)
sfExportAll()
prob <- sfSapply( 1:(n-1),
function(i){
ensemble(starts[i,], initial=ends[i,], pars = pars, dt = dt, seed =seed[i], reps = reps)$probs[5]
})
prob
# -sum(log(prob))
}