-
Notifications
You must be signed in to change notification settings - Fork 1
/
utils.R
47 lines (39 loc) · 1.11 KB
/
utils.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
# Function that does MCMC
do_mcmc = function(particle, log_likelihood, threshold)
{
accepted = 0
for(j in 1:mcmc_steps)
{
# Generate proposal by perturbing one or several of the underlying
# coordinates
new = particle
num = 1
if(runif(1) < 0.5)
num = floor(num_params**runif(1))
for(k in 1:num)
{
which = sample(1:num_params, 1)
scale = 10.0**(1.0 - 3*abs(rt(1, df=2)))
new[which] = new[which] + rnorm(1, sd=scale)
new[which] = new[which] %% 1.0
}
# Accept the proposal if the new particle is above the threshold
logL_new = log_likelihood(us_to_params(new))
if(logL_new >= threshold)
{
particle = new
log_likelihood = logL_new
accepted = accepted + 1
}
}
return(list(particle=particle,
log_likelihood=log_likelihood, accepted=accepted))
}
# Useful function
logsumexp = function(xs)
{
biggest = max(xs)
xs = xs - biggest
result = log(sum(exp(xs))) + biggest
return(result)
}