@@ -27,7 +27,7 @@ Holds state of algorithm
27
27
"""
28
28
mutable struct SABCstate
29
29
ϵ:: Vector{Float64} # epsilon = temperature
30
- algorithm:: Symbol # the algorithm used
30
+ algorithm:: Symbol # the algorithm used
31
31
32
32
# Containers to store trajectories
33
33
ϵ_history:: Vector{Vector{Float64}}
@@ -50,6 +50,8 @@ Holds results from a SABC run with fields:
50
50
- `state`: state of algorithm
51
51
52
52
The history of ϵ can be accessed with the field `state.ϵ_history`.
53
+ The history of ρ can be accessed with the field `state.ρ_history`.
54
+ The history of u can be accessed with the field `state.u_history`.
53
55
"""
54
56
struct SABCresult{T, S}
55
57
population:: Vector{T}
@@ -66,17 +68,18 @@ function show(io::Base.IO, s::SABCresult)
66
68
mean_u = round (mean (s. u), sigdigits = 4 )
67
69
acc_rate = round (s. state. n_accept / (s. state. n_simulation - n_particles),
68
70
sigdigits = 4 )
69
-
70
71
println (io, " Approximate posterior sample with $n_particles particles:" )
71
72
println (io, " - algorithm: :$(s. state. algorithm) " )
72
73
println (io, " - simulations used: $(s. state. n_simulation) " )
73
74
println (io, " - number of population updates: $(s. state. n_population_updates) " )
74
75
println (io, " - average transformed distance: $mean_u " )
75
76
println (io, " - ϵ: $(round .(s. state. ϵ, sigdigits= 4 )) " )
76
- println (io, " - population resampling : $(s. state. n_resampling) " )
77
+ println (io, " - number of population resamplings : $(s. state. n_resampling) " )
77
78
println (io, " - acceptance rate: $(acc_rate) " )
78
79
println (io, " The sample can be accessed with the field `population`." )
79
80
println (io, " The history of ϵ can be accessed with the field `state.ϵ_history`." )
81
+ println (io, " The history of ρ can be accessed with the field `state.ρ_history`." )
82
+ println (io, " The history of u can be accessed with the field `state.u_history`." )
80
83
end
81
84
82
85
# -------------------------------------------
@@ -93,7 +96,7 @@ function update_epsilon_single_stat(ū, v)
93
96
end
94
97
95
98
"""
96
- Update ϵ for multi-epsilons. See eq(19-20) in Albert et al. (in preparatory).
99
+ Update ϵ for multi-epsilons. See eq(19-20) in Albert et al. (in preparation)
97
100
"""
98
101
function update_epsilon_multi_stats (u, v)
99
102
n = size (u, 2 ) # number of statistics
120
123
Resample population
121
124
"""
122
125
function resample_population (population, u, δ)
123
-
124
126
n = length (population)
125
- u_means = mean (u, dims= 1 )
126
- w = exp .(- sum (u[:,i] .* δ ./ u_means [i] for i in 1 : size (u, 2 )))
127
+ ū = mean (u, dims= 1 )
128
+ w = exp .(- sum (u[:,i] .* δ ./ ū [i] for i in 1 : size (u, 2 )))
127
129
# Choose indexes based on weights w
128
130
idx_resampled = sample (1 : n, weights (w), n, replace= true )
129
131
# Apply selected indexes to population and u's
@@ -160,7 +162,7 @@ proposal(θ, Σ::Float64) = θ + rand(Normal(0, sqrt(Σ)))
160
162
## Arguments
161
163
See docs for `sabc`
162
164
163
- ## Return
165
+ ## Returns
164
166
- An object of type `SABCresult`
165
167
"""
166
168
function initialization (f_dist, prior:: Distribution , args... ;
@@ -194,7 +196,6 @@ function initialization(f_dist, prior::Distribution, args...;
194
196
end
195
197
ρ_history = [[mean (ic) for ic in eachcol (distances_prior)]]
196
198
197
-
198
199
# ------------------
199
200
# Estimate the cdf of ρ under the prior
200
201
@@ -223,10 +224,12 @@ function initialization(f_dist, prior::Distribution, args...;
223
224
u_history = [[mean (ic) for ic in eachcol (u)]]
224
225
ϵ_history = [copy (ϵ)]
225
226
227
+ # ------------------
228
+ # estimate jump covariance
226
229
Σ_jump = estimate_jump_covariance (population, β)
227
230
228
231
# ---------------------
229
- # Collect parameters and states of the algorithm
232
+ # Collect parameters and state of the algorithm
230
233
231
234
n_simulation = n_particles # N.B.: we consider only n_particles draws from the prior
232
235
# and neglect the first call to f_dist ('initialization of containers')
@@ -290,8 +293,8 @@ function update_population!(population_state::SABCresult, f_dist, prior, args...
290
293
291
294
n_particles = length (population)
292
295
293
- n_population_updates = n_simulation ÷ n_particles # populations update = total simulations / particles
294
- n_updates = n_population_updates * n_particles # number of calls to `f_dist`
296
+ n_population_updates = n_simulation ÷ n_particles # population updates = total simulations / particles
297
+ n_updates = n_population_updates * n_particles # number of calls to `f_dist` (individual particle updates)
295
298
last_checkpoint_epsilon = 0 # set checkpoint counter to zero
296
299
297
300
# to estimate ETA
@@ -303,6 +306,7 @@ function update_population!(population_state::SABCresult, f_dist, prior, args...
303
306
pmeter = Progress (n_population_updates; desc = " $n_population_updates population updates:" ,
304
307
output = stderr , enabled = show_progressbar)
305
308
generate_showvalues (ϵ, u) = () -> [(" ϵ" , round .(ϵ, sigdigits= 4 )), (" average transformed distance" , round .(mean (u), sigdigits= 4 ))]
309
+
306
310
for ix in 1 : n_population_updates
307
311
308
312
# ----------------------------------------------------------
0 commit comments