From a373707a5e1c3a6c2ca86c54164c1d59686b8e62 Mon Sep 17 00:00:00 2001 From: Simone Ulzega Date: Thu, 26 Sep 2024 11:42:46 +0200 Subject: [PATCH] minor things, typos and alike --- src/SimulatedAnnealingABC.jl | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/src/SimulatedAnnealingABC.jl b/src/SimulatedAnnealingABC.jl index c3fedce..d056eea 100644 --- a/src/SimulatedAnnealingABC.jl +++ b/src/SimulatedAnnealingABC.jl @@ -27,7 +27,7 @@ Holds state of algorithm """ mutable struct SABCstate ϵ::Vector{Float64} # epsilon = temperature - algorithm::Symbol # the algorithm used + algorithm::Symbol # the algorithm used # Containers to store trajectories ϵ_history::Vector{Vector{Float64}} @@ -50,6 +50,8 @@ Holds results from a SABC run with fields: - `state`: state of algorithm The history of ϵ can be accessed with the field `state.ϵ_history`. +The history of ρ can be accessed with the field `state.ρ_history`. +The history of u can be accessed with the field `state.u_history`. """ struct SABCresult{T, S} population::Vector{T} @@ -66,17 +68,18 @@ function show(io::Base.IO, s::SABCresult) mean_u = round(mean(s.u), sigdigits = 4) acc_rate = round(s.state.n_accept / (s.state.n_simulation - n_particles), sigdigits = 4) - println(io, "Approximate posterior sample with $n_particles particles:") println(io, " - algorithm: :$(s.state.algorithm)") println(io, " - simulations used: $(s.state.n_simulation)") println(io, " - number of population updates: $(s.state.n_population_updates)") println(io, " - average transformed distance: $mean_u") println(io, " - ϵ: $(round.(s.state.ϵ, sigdigits=4))") - println(io, " - population resampling: $(s.state.n_resampling)") + println(io, " - number of population resamplings: $(s.state.n_resampling)") println(io, " - acceptance rate: $(acc_rate)") println(io, "The sample can be accessed with the field `population`.") println(io, "The history of ϵ can be accessed with the field `state.ϵ_history`.") + println(io, "The history of ρ can be accessed with the field `state.ρ_history`.") + println(io, "The history of u can be accessed with the field `state.u_history`.") end # ------------------------------------------- @@ -93,7 +96,7 @@ function update_epsilon_single_stat(ū, v) end """ -Update ϵ for multi-epsilons. See eq(19-20) in Albert et al. (in preparatory). +Update ϵ for multi-epsilons. See eq(19-20) in Albert et al. (in preparation) """ function update_epsilon_multi_stats(u, v) n = size(u, 2) # number of statistics @@ -120,10 +123,9 @@ end Resample population """ function resample_population(population, u, δ) - n = length(population) - u_means = mean(u, dims=1) - w = exp.(-sum(u[:,i] .* δ ./ u_means[i] for i in 1:size(u, 2))) + ū = mean(u, dims=1) + w = exp.(-sum(u[:,i] .* δ ./ ū[i] for i in 1:size(u, 2))) # Choose indexes based on weights w idx_resampled = sample(1:n, weights(w), n, replace=true) # Apply selected indexes to population and u's @@ -160,7 +162,7 @@ proposal(θ, Σ::Float64) = θ + rand(Normal(0, sqrt(Σ))) ## Arguments See docs for `sabc` -## Return +## Returns - An object of type `SABCresult` """ function initialization(f_dist, prior::Distribution, args...; @@ -194,7 +196,6 @@ function initialization(f_dist, prior::Distribution, args...; end ρ_history = [[mean(ic) for ic in eachcol(distances_prior)]] - # ------------------ # Estimate the cdf of ρ under the prior @@ -223,10 +224,12 @@ function initialization(f_dist, prior::Distribution, args...; u_history = [[mean(ic) for ic in eachcol(u)]] ϵ_history = [copy(ϵ)] + # ------------------ + # estimate jump covariance Σ_jump = estimate_jump_covariance(population, β) # --------------------- - # Collect parameters and states of the algorithm + # Collect parameters and state of the algorithm n_simulation = n_particles # N.B.: we consider only n_particles draws from the prior # 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... n_particles = length(population) - n_population_updates = n_simulation ÷ n_particles # populations update = total simulations / particles - n_updates = n_population_updates * n_particles # number of calls to `f_dist` + n_population_updates = n_simulation ÷ n_particles # population updates = total simulations / particles + n_updates = n_population_updates * n_particles # number of calls to `f_dist` (individual particle updates) last_checkpoint_epsilon = 0 # set checkpoint counter to zero # to estimate ETA @@ -303,6 +306,7 @@ function update_population!(population_state::SABCresult, f_dist, prior, args... pmeter = Progress(n_population_updates; desc = "$n_population_updates population updates:", output = stderr, enabled = show_progressbar) generate_showvalues(ϵ, u) = () -> [("ϵ", round.(ϵ, sigdigits=4)), ("average transformed distance", round.(mean(u), sigdigits=4))] + for ix in 1:n_population_updates # ----------------------------------------------------------