/
metropolis_tuned.py
75 lines (55 loc) · 2.27 KB
/
metropolis_tuned.py
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
def metropolis_tuned(n_iterations, initial_values, prop_var=1,
tune_for=None, tune_interval=100):
n_params = len(initial_values)
# Initial proposal standard deviations
prop_sd = [prop_var] * n_params
# Initialize trace for parameters
trace = np.empty((n_iterations+1, n_params))
# Set initial values
trace[0] = initial_values
# Initialize acceptance counts
accepted = [0]*n_params
# Calculate joint posterior for initial values
current_log_prob = calc_posterior(*trace[0])
if tune_for is None:
tune_for = n_iterations/2
for i in range(n_iterations):
if not i%1000: print('Iteration %i' % i)
# Grab current parameter values
current_params = trace[i]
for j in range(n_params):
# Get current value for parameter j
p = trace[i].copy()
# Propose new value
if j==2:
# Ensure tau is positive
theta = np.exp(rnorm(np.log(current_params[j]), prop_sd[j]))
else:
theta = rnorm(current_params[j], prop_sd[j])
# Insert new value
p[j] = theta
# Calculate log posterior with proposed value
proposed_log_prob = calc_posterior(*p)
# Log-acceptance rate
alpha = proposed_log_prob - current_log_prob
# Sample a uniform random variate
u = runif()
# Test proposed value
if np.log(u) < alpha:
# Accept
trace[i+1,j] = theta
current_log_prob = proposed_log_prob
accepted[j] += 1
else:
# Reject
trace[i+1,j] = trace[i,j]
# Tune every 100 iterations
if (not (i+1) % tune_interval) and (i < tune_for):
# Calculate aceptance rate
acceptance_rate = (1.*accepted[j])/tune_interval
if acceptance_rate<0.2:
prop_sd[j] *= 0.9
elif acceptance_rate>0.5:
prop_sd[j] *= 1.1
accepted[j] = 0
return trace[tune_for:], accepted