From f36753c68a27d9270c8c760e53b931f11692f36c Mon Sep 17 00:00:00 2001 From: William Chen Date: Thu, 18 Mar 2021 13:05:07 -0400 Subject: [PATCH] Resolve logpost logprior confusion (column of cloud.particles is logprior, not logpost) --- src/initialization.jl | 32 ++++++++++++++++---------------- src/mutation.jl | 6 +++--- src/particle.jl | 40 +++++++++++++++++++--------------------- src/smc_main.jl | 6 +++--- 4 files changed, 41 insertions(+), 43 deletions(-) diff --git a/src/initialization.jl b/src/initialization.jl index c1f15db..b4f919d 100644 --- a/src/initialization.jl +++ b/src/initialization.jl @@ -6,7 +6,7 @@ one_draw(loglikelihood::Function, parameters::ParameterVector{U}, ``` Finds and returns one valid draw from parameter distribution, along with its -log likelihood and log posterior. +log likelihood and log prior. Set `regime_switching` to true if there are regime-switching parameters. Otherwise, not all the values of the @@ -26,7 +26,7 @@ function one_draw(loglikelihood::Function, parameters::ParameterVector{U}, success = false draw = vec(rand(parameters, 1, regime_switching = regime_switching, toggle = toggle)) - draw_loglh = draw_logpost = 0.0 + draw_loglh = draw_logprior = 0.0 while !success try @@ -38,16 +38,16 @@ function one_draw(loglikelihood::Function, parameters::ParameterVector{U}, toggle_regime!(parameters, 1) end - draw_logpost = prior(parameters) + draw_logprior = prior(parameters) if (draw_loglh == -Inf) || (draw_loglh === NaN) - draw_loglh = draw_logpost = -Inf + draw_loglh = draw_logprior = -Inf end catch err if isa(err, ParamBoundsError) || isa(err, SingularException) || isa(err, LinearAlgebra.LAPACKException) || isa(err, PosDefException) || isa(err, DomainError) - draw_loglh = draw_logpost = -Inf + draw_loglh = draw_logprior = -Inf else throw(err) end @@ -59,7 +59,7 @@ function one_draw(loglikelihood::Function, parameters::ParameterVector{U}, success = true end end - return vector_reshape(draw, draw_loglh, draw_logpost) + return vector_reshape(draw, draw_loglh, draw_logprior) end """ @@ -70,7 +70,7 @@ function initial_draw!(loglikelihood::Function, parameters::ParameterVector{U}, ``` Draw from a general starting distribution (set by default to be from the prior) to -initialize the SMC algorithm. Returns a tuple (logpost, loglh) and modifies the +initialize the SMC algorithm. Returns a tuple (loglh, logprior) and modifies the particle objects in the particle cloud in place. Set `regime_switching` to true if @@ -99,8 +99,8 @@ function initial_draw!(loglikelihood::Function, parameters::ParameterVector{U}, @everywhere one_draw_closure() = one_draw(loglikelihood, parameters, data, regime_switching = regime_switching, toggle = toggle) # ========================================================================= - # For each particle, finds valid parameter draw and returns loglikelihood & posterior - draws, loglh, logpost = if parallel + # For each particle, finds valid parameter draw and returns loglikelihood & prior + draws, loglh, logprior = if parallel @sync @distributed (vector_reduce) for i in 1:n_parts one_draw_closure() end @@ -110,7 +110,7 @@ function initial_draw!(loglikelihood::Function, parameters::ParameterVector{U}, update_draws!(c, draws) update_loglh!(c, vec(loglh)) - update_logpost!(c, vec(logpost)) + update_logprior!(c, vec(logprior)) update_old_loglh!(c, zeros(n_parts)) # Need to call `set_weights` as opposed to `update_weights` @@ -124,7 +124,7 @@ draw_likelihood(loglikelihood::Function, parameters::ParameterVector{U}, data::Matrix{Float64}, draw::Vector{Float64}; toggle::Bool = true) where {U<:Number} ``` -Computes likelihood of a particular parameter draw; returns loglh and logpost. +Computes likelihood of a particular parameter draw; returns loglh and logprior. """ function draw_likelihood(loglikelihood::Function, parameters::ParameterVector{U}, data::Matrix{Float64}, draw::Vector{Float64}; @@ -134,8 +134,8 @@ function draw_likelihood(loglikelihood::Function, parameters::ParameterVector{U} if toggle ModelConstructors.toggle_regime!(parameters, 1) end - logpost = prior(parameters) - return scalar_reshape(loglh, logpost) + logprior = prior(parameters) + return scalar_reshape(loglh, logprior) end """ @@ -147,7 +147,7 @@ initialize_likelihoods!(loglikelihood::Function, parameters::ParameterVector{U}, ``` This function is made for transfering the log-likelihood values saved in the Cloud from a previous estimation to each particle's respective old_loglh -field, and for evaluating/saving the likelihood and posterior at the new data, which +field, and for evaluating/saving the likelihood and prior at the new data, which here is just the argument, data. """ function initialize_likelihoods!(loglikelihood::Function, parameters::ParameterVector{U}, @@ -174,7 +174,7 @@ function initialize_likelihoods!(loglikelihood::Function, parameters::ParameterV # TODO: handle when the likelihood with new data cannot be evaluated (returns -Inf), # even if the likelihood was not -Inf prior to incorporating new data - loglh, logpost = if parallel + loglh, logprior = if parallel @sync @distributed (scalar_reduce) for i in 1:n_parts draw_likelihood_closure(draws[i, :]) end @@ -182,7 +182,7 @@ function initialize_likelihoods!(loglikelihood::Function, parameters::ParameterV scalar_reduce([draw_likelihood_closure(draws[i, :]) for i in 1:n_parts]...) end update_loglh!(c, loglh) - update_logpost!(c, logpost) + update_logprior!(c, logprior) end """ diff --git a/src/mutation.jl b/src/mutation.jl index 98caa8d..4148999 100755 --- a/src/mutation.jl +++ b/src/mutation.jl @@ -32,7 +32,7 @@ Execute one proposed move of the Metropolis-Hastings algorithm for a given param - `α::Float64`: The mixing proportion - `n_mh_steps::Int`: Number of Metropolis Hastings steps to attempt - `old_data::Matrix{Float64}`: The matrix of old data to be used in calculating the - old_loglh, old_logpost in time tempering + old_loglh, old_logprior in time tempering - `regime_switching::Bool`: Set to true if there are regime-switching parameters. Otherwise, not all the values of the regimes will be used or saved. @@ -83,7 +83,7 @@ function mutation(loglikelihood::Function, parameters::ParameterVector{U}, prior_new = like_new = like_old_data = -Inf try update!(parameters, para_new) - para_new = [θ.value for θ in parameters] + #=para_new = [θ.value for θ in parameters] i = 0 for para in parameters if !isempty(para.regimes) @@ -94,7 +94,7 @@ function mutation(loglikelihood::Function, parameters::ParameterVector{U}, end end end - end + end=# prior_new = prior(parameters) like_new = loglikelihood(parameters, data) #, regime_switching = regime_switching) diff --git a/src/particle.jl b/src/particle.jl index e7290d1..bf7c0fd 100644 --- a/src/particle.jl +++ b/src/particle.jl @@ -53,8 +53,7 @@ Find correct indices for accessing columns of cloud array. """ ind_para_end(N::Int) = N-5 ind_loglh(N::Int) = N-4 -ind_logpost(N::Int) = N-3 -ind_logprior(N::Int) = N-3 # TODO: Fix logprior/logpost shenanigans +ind_logprior(N::Int) = N-3 ind_old_loglh(N::Int) = N-2 ind_accept(N::Int) = N-1 ind_weight(N::Int) = N @@ -166,7 +165,6 @@ function get_logprior(c::Cloud) Returns Vector{Float64}(n_parts) of log-prior of particles in cloud. """ function get_logprior(c::Matrix{Float64}) - #TODO: Fix logpost/logprior confusion return c[:, ind_logprior(size(c,2))] end function get_logprior(c::Cloud) @@ -278,20 +276,20 @@ end """ ``` -function update_logpost!(c::Matrix{Float64}, incweight::Vector{Float64}) -function update_logpost!(c::Cloud, logpost::Vector{Float64}) +function update_logprior!(c::Matrix{Float64}, incweight::Vector{Float64}) +function update_logprior!(c::Cloud, logprior::Vector{Float64}) ``` -Update log-posterior in cloud. +Update log-prior in cloud. """ -function update_logpost!(c::Matrix{Float64}, logpost::Vector{Float64}) - @assert size(c, 1) == length(logpost) "Dimensional mismatch" - N = ind_logpost(size(c,2)) - for i=1:length(logpost) - c[i, N] = logpost[i] +function update_logprior!(c::Matrix{Float64}, logprior::Vector{Float64}) + @assert size(c, 1) == length(logprior) "Dimensional mismatch" + N = ind_logprior(size(c,2)) + for i=1:length(logprior) + c[i, N] = logprior[i] end end -function update_logpost!(c::Cloud, logpost::Vector{Float64}) - update_logpost!(c.particles, logpost) +function update_logprior!(c::Cloud, logprior::Vector{Float64}) + update_logprior!(c.particles, logprior) end @@ -347,18 +345,18 @@ end """ ``` function update_mutation!(p::Vector{Float64}, para::Vector{Float64}, - like::Float64, post::Float64, old_like::Float64, + like::Float64, prior::Float64, old_like::Float64, accept::Float64) ``` -Update a particle's parameter vector, log-likelihood, log-posteriod, +Update a particle's parameter vector, log-likelihood, log-prior, old log-likelihood, and acceptance rate at the end of mutation. """ function update_mutation!(p::Vector{Float64}, para::Vector{Float64}, like::Float64, - post::Float64, old_like::Float64, accept::Float64) + prior::Float64, old_like::Float64, accept::Float64) N = length(p) p[1:ind_para_end(N)] = para p[ind_loglh(N)] = like - p[ind_logpost(N)] = post + p[ind_logprior(N)] = prior p[ind_old_loglh(N)] = old_like p[ind_accept(N)] = accept end @@ -496,7 +494,7 @@ function split_cloud(filename::String, n_pieces::Int) Ws[i] = W[inds, :] clouds[i] = Cloud(n_para, npart_small) - # Since loglh, logpost, old_loglh are all stored in cloud.particles, this updates them all too! + # Since loglh, logprior, old_loglh are all stored in cloud.particles, this updates them all too! SMC.update_cloud!(clouds[i], cloud.particles[inds, :]) clouds[i].ESS = cloud.ESS @@ -537,7 +535,7 @@ function join_cloud(filename::String, n_pieces::Int) particles = Matrix{Float64}(undef, 0, n_para) loglhs = Vector{Float64}(undef, 0) old_loglhs = Vector{Float64}(undef, 0) - logposts = Vector{Float64}(undef, 0) + logpriors = Vector{Float64}(undef, 0) w = Matrix{Float64}(undef, 0, n_stage) W = Matrix{Float64}(undef, 0, n_stage) @@ -545,12 +543,12 @@ function join_cloud(filename::String, n_pieces::Int) particles = vcat(particles, clouds[i].particles) loglhs = vcat(loglhs, SMC.get_loglh(clouds[i])) old_loglhs = vcat(old_loglhs, SMC.get_old_loglh(clouds[i])) - logposts = vcat(logposts, SMC.get_logpost(clouds[i])) + logpriors = vcat(logpriors, SMC.get_logprior(clouds[i])) w = vcat(w, ws[i]) W = vcat(W, Ws[i]) end - # Since loglh, logpost, old_loglh are all stored in cloud.particles, this updates them all too! + # Since loglh, logprior, old_loglh are all stored in cloud.particles, this updates them all too! SMC.update_cloud!(cloud, particles) cloud.ESS = clouds[1].ESS diff --git a/src/smc_main.jl b/src/smc_main.jl index 3dd3dd2..75ced68 100644 --- a/src/smc_main.jl +++ b/src/smc_main.jl @@ -255,7 +255,7 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr bridge_cloud = Cloud(n_para, n_to_resample) update_cloud!(bridge_cloud, cloud.particles[new_inds, :]) update_loglh!(bridge_cloud, get_loglh(cloud)[new_inds]) - update_logpost!(bridge_cloud, get_logpost(cloud)[new_inds]) + update_logprior!(bridge_cloud, get_logprior(cloud)[new_inds]) update_old_loglh!(bridge_cloud, get_old_loglh(cloud)[new_inds]) # Make a cloud by drawing from the prior @@ -267,7 +267,7 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr cloud = Cloud(n_para, n_to_resample + n_from_prior) update_cloud!(cloud, vcat(bridge_cloud.particles, prior_cloud.particles)) update_loglh!(cloud, vcat(get_loglh(bridge_cloud), get_loglh(prior_cloud))) - update_logpost!(cloud, vcat(get_logpost(bridge_cloud), get_logpost(prior_cloud))) + update_logprior!(cloud, vcat(get_logprior(bridge_cloud), get_logprior(prior_cloud))) update_old_loglh!(cloud, vcat(get_old_loglh(bridge_cloud), get_old_loglh(prior_cloud))) else cloud = bridge_cloud @@ -286,7 +286,7 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr elseif continue_intermediate cloud = load(loadpath, "cloud") else - # Instantiating Cloud object, update draws, loglh, & logpost + # Instantiating Cloud object, update draws, loglh, & logprior initial_draw!(loglikelihood, parameters, data, cloud; parallel = parallel, regime_switching = regime_switching, toggle = toggle) initialize_cloud_settings!(cloud; tempered_update = tempered_update,