Resolve logpost logprior confusion (column of cloud.particles is logp…

…rior, not logpost)
FRBNY-DSGE · Mar 18, 2021 · f36753c · f36753c
1 parent a6ade4f
commit f36753c
Show file tree

Hide file tree

Showing 4 changed files with 41 additions and 43 deletions.
diff --git 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,15 +174,15 @@ 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
     else
         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
@@ -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
@@ -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,20 +535,20 @@ 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)
 
     for i = 1:n_pieces
         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
@@ -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,