Skip to content

Commit

Permalink
Fix likelihood for linear reg model. Update reference output. Add che…
Browse files Browse the repository at this point in the history 
…ck that parameter estimates are close to true values.
  • Loading branch information
@chenwilliam77
chenwilliam77 committed Dec 18, 2020
1 parent 7bb05ee commit b5ed5a4
Show file tree
Hide file tree
Showing 9 changed files with 220 additions and 81 deletions.
68 changes: 49 additions & 19 deletions src/initialization.jl
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,8 @@
"""
```
function one_draw(loglikelihood::Function, parameters::ParameterVector{U},
data::Matrix{Float64}; regime_switching::Bool = false) where {U<:Number}
one_draw(loglikelihood::Function, parameters::ParameterVector{U},
data::Matrix{Float64}; regime_switching::Bool = false,
toggle::Bool = true) where {U<:Number}
```
Finds and returns one valid draw from parameter distribution, along with its
Expand All @@ -10,19 +11,33 @@ log likelihood and log posterior.
Set `regime_switching` to true if
there are regime-switching parameters. Otherwise, not all the values of the
regimes will be used or saved.
Set `toggle` to false if, after calculating the loglikelihood,
the values in the fields of every parameter in `parameters`
are set to their regime 1 values. The regime-switching version of `rand`
requires that the fields of all parameters take their regime 1 values,
or else sampling may be wrong. The default is `true` as a safety, but
if speed is a paramount concern, setting `toggle = true` will avoid
unnecessary computations.
"""
function one_draw(loglikelihood::Function, parameters::ParameterVector{U},
data::Matrix{Float64}; regime_switching::Bool = false) where {U<:Number}
data::Matrix{Float64}; regime_switching::Bool = false,
toggle::Bool = true) where {U<:Number}
success = false
draw = vec(rand(parameters, 1, regime_switching = regime_switching))
draw = vec(rand(parameters, 1, regime_switching = regime_switching, toggle = toggle))

draw_loglh = draw_logpost = 0.0

while !success
try
update!(parameters, draw)

draw_loglh = loglikelihood(parameters, data)
draw_loglh = loglikelihood(parameters, data)

if toggle
toggle_regime!(parameters, 1)
end

draw_logpost = prior(parameters)

if (draw_loglh == -Inf) || (draw_loglh === NaN)
Expand All @@ -39,7 +54,7 @@ function one_draw(loglikelihood::Function, parameters::ParameterVector{U},
end

if any(isinf.(draw_loglh))
draw = vec(rand(parameters, 1, regime_switching = regime_switching))
draw = vec(rand(parameters, 1, regime_switching = regime_switching, toggle = false))
else
success = true
end
Expand All @@ -51,7 +66,7 @@ end
```
function initial_draw!(loglikelihood::Function, parameters::ParameterVector{U},
data::Matrix{Float64}, c::Cloud; parallel::Bool = false,
regime_switching::Bool = false) where {U<:Number}
regime_switching::Bool = false, toggle::Bool = true) where {U<:Number}
```
Draw from a general starting distribution (set by default to be from the prior) to
Expand All @@ -61,19 +76,27 @@ particle objects in the particle cloud in place.
Set `regime_switching` to true if
there are regime-switching parameters. Otherwise, not all the values of the
regimes will be used or saved.
Set `toggle` to false if, after calculating the loglikelihood,
the values in the fields of every parameter in `parameters`
are set to their regime 1 values. The regime-switching version of `rand`
requires that the fields of all parameters take their regime 1 values,
or else sampling may be wrong. The default is `true` as a safety, but
if speed is a paramount concern, setting `toggle = true` will avoid
unnecessary computations.
"""
function initial_draw!(loglikelihood::Function, parameters::ParameterVector{U},
data::Matrix{Float64}, c::Cloud; parallel::Bool = false,
regime_switching::Bool = false) where {U<:Number}
regime_switching::Bool = false, toggle::Bool = true) where {U<:Number}
n_parts = length(c)

# ================== Define closure on one_draw function ==================
sendto(workers(), loglikelihood = loglikelihood)
sendto(workers(), parameters = parameters)
sendto(workers(), data = data)

one_draw_closure() = one_draw(loglikelihood, parameters, data, regime_switching = regime_switching)
@everywhere one_draw_closure() = one_draw(loglikelihood, parameters, data, regime_switching = regime_switching)
one_draw_closure() = one_draw(loglikelihood, parameters, data, regime_switching = regime_switching, toggle = toggle)
@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
Expand All @@ -97,24 +120,30 @@ end

"""
```
function draw_likelihood(loglikelihood::Function, parameters::ParameterVector{U},
data::Matrix{Float64}, draw::Vector{Float64}) where {U<:Number}
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.
"""
function draw_likelihood(loglikelihood::Function, parameters::ParameterVector{U},
data::Matrix{Float64}, draw::Vector{Float64}) where {U<:Number}
data::Matrix{Float64}, draw::Vector{Float64};
toggle::Bool = true) where {U<:Number}
update!(parameters, draw)
loglh = loglikelihood(parameters, data)
if toggle
ModelConstructors.toggle_regime!(parameters, 1)
end
logpost = prior(parameters)
return scalar_reshape(loglh, logpost)
end

"""
```
function initialize_likelihoods!(loglikelihood::Function, parameters::ParameterVector{U},
data::Matrix{Float64}, c::Cloud;
parallel::Bool = false) where {U<:Number}
initialize_likelihoods!(loglikelihood::Function, parameters::ParameterVector{U},
data::Matrix{Float64}, c::Cloud;
parallel::Bool = false,
toggle::Bool = true) where {U<:Number}
```
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
Expand All @@ -123,7 +152,8 @@ here is just the argument, data.
"""
function initialize_likelihoods!(loglikelihood::Function, parameters::ParameterVector{U},
data::Matrix{Float64}, c::Cloud;
parallel::Bool = false) where {U<:Number}
parallel::Bool = false,
toggle::Bool = true) where {U<:Number}
n_parts = length(c)
draws = get_vals(c; transpose = false)

Expand All @@ -136,10 +166,10 @@ function initialize_likelihoods!(loglikelihood::Function, parameters::ParameterV
sendto(workers(), data = data)

draw_likelihood_closure(draw::Vector{Float64}) = draw_likelihood(loglikelihood, parameters,
data, draw)
data, draw, toggle = toggle)
@everywhere draw_likelihood_closure(draw::Vector{Float64}) = draw_likelihood(loglikelihood,
parameters,
data, draw)
data, draw, toggle = toggle)
# =========================================================================

# TODO: handle when the likelihood with new data cannot be evaluated (returns -Inf),
Expand Down
15 changes: 14 additions & 1 deletion src/mutation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,10 @@ Execute one proposed move of the Metropolis-Hastings algorithm for a given param
- `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.
- `toggle::Bool = true`: Flag for resetting the fields of parameter values to regime 1 anytime
the loglikelihood is computed. The regime-switching version of SMC assumes at various points
that this resetting occurs. If speed is important, then ensure that the fields of parameters
take their regime 1 values at the end of the loglikelihood computation and set `toggle = false`.
### Outputs:
- `p::Vector{Float64}`: An updated particle containing updated parameter values,
Expand All @@ -48,7 +52,8 @@ function mutation(loglikelihood::Function, parameters::ParameterVector{U},
blocks_free::Vector{Vector{Int}}, blocks_all::Vector{Vector{Int}},
ϕ_n::S, ϕ_n1::S; c::S = 1., α::S = 1., n_mh_steps::Int = 1,
old_data::T = T(undef, size(data, 1), 0),
regime_switching::Bool = false) where {S<:AbstractFloat,T<:AbstractMatrix, U<:Number}
regime_switching::Bool = false,
toggle::Bool = true) where {S<:AbstractFloat,T<:AbstractMatrix, U<:Number}

step_prob = rand() # Draw initial step probability

Expand Down Expand Up @@ -94,12 +99,20 @@ function mutation(loglikelihood::Function, parameters::ParameterVector{U},
prior_new = prior(parameters)
like_new = loglikelihood(parameters, data) #, regime_switching = regime_switching)

if toggle
toggle_regime!(parameters, 1)
end

if like_new == -Inf
prior_new = like_old_data = -Inf
end

like_old_data = isempty(old_data) ? 0. : loglikelihood(parameters, old_data)

if toggle && isempty(old_data)
toggle_regime!(parameters, 1)
end

catch err
if isa(err, ParamBoundsError) || isa(err, LinearAlgebra.LAPACKException) ||
isa(err, PosDefException) || isa(err, SingularException) ||
Expand Down
44 changes: 28 additions & 16 deletions src/smc_main.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,10 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr
the Hessian is negative definite).
- `regime_switching::Bool = false`: Flag if there are regime-switching parameters. Otherwise, not all the values of the
regimes will be used or saved.
- `toggle::Bool = true`: Flag for resetting the fields of parameter values to regime 1 anytime
the loglikelihood is computed. The regime-switching version of SMC assumes at various points
that this resetting occurs. If speed is important, then ensure that the fields of parameters
take their regime 1 values at the end of the loglikelihood computation and set `toggle = false`.
### Outputs
Expand Down Expand Up @@ -141,7 +145,8 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr
continue_intermediate::Bool = false,
intermediate_stage_start::Int = 0,
tempered_update_prior_weight::S = 0.0,
regime_switching::Bool = false) where {S<:AbstractFloat, U<:Number}
regime_switching::Bool = false,
toggle::Bool = true) where {S<:AbstractFloat, U<:Number}

########################################################################################
### Settings
Expand All @@ -157,16 +162,19 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr
blocks_free::Vector{Vector{Int64}}, blocks_all::Vector{Vector{Int64}},
ϕ_n::S, ϕ_n1::S; c::S = 1.0, α::S = 1.0, n_mh_steps::Int = 1,
old_data::T = Matrix{S}(undef, size(data, 1), 0),
regime_switching::Bool = false) where {S<:Float64, T<:Matrix}
regime_switching::Bool = false, toggle::Bool = true) where {S<:Float64, T<:Matrix}
return mutation(loglikelihood, parameters, data, p, d_μ, d_Σ, n_free_para, blocks_free, blocks_all,
ϕ_n, ϕ_n1; c = c, α = α, n_mh_steps = n_mh_steps, old_data = old_data, regime_switching = regime_switching)
ϕ_n, ϕ_n1; c = c, α = α, n_mh_steps = n_mh_steps, old_data = old_data, regime_switching = regime_switching,
toggle = toggle)
end
@everywhere function mutation_closure(p::Vector{S}, d_μ::Vector{S}, d_Σ::Matrix{S},
blocks_free::Vector{Vector{Int64}}, blocks_all::Vector{Vector{Int64}}, n_free_para::Int,
ϕ_n::S, ϕ_n1::S; c::S = 1.0, α::S = 1.0, n_mh_steps::Int = 1,
old_data::T = Matrix{S}(undef, size(data, 1), 0), regime_switching::Bool = false) where {S<:Float64, T<:Matrix}
blocks_free::Vector{Vector{Int64}}, blocks_all::Vector{Vector{Int64}}, n_free_para::Int,
ϕ_n::S, ϕ_n1::S; c::S = 1.0, α::S = 1.0, n_mh_steps::Int = 1,
old_data::T = Matrix{S}(undef, size(data, 1), 0), regime_switching::Bool = false,
toggle::Bool = true) where {S<:Float64, T<:Matrix}
return mutation(loglikelihood, parameters, data, p, d_μ, d_Σ, blocks_free, blocks_all, n_free_para,
ϕ_n, ϕ_n1; c = c, α = α, n_mh_steps = n_mh_steps, old_data = old_data, regime_switching = regime_switching)
ϕ_n, ϕ_n1; c = c, α = α, n_mh_steps = n_mh_steps, old_data = old_data,
regime_switching = regime_switching, toggle = toggle)
end

# Check that if there's a tempered update, old and current vintages are different
Expand All @@ -189,9 +197,7 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr
for para in parameters
if !isempty(para.regimes)
for (ind, val) in para.regimes[:value]
if ind == 1
nothing
else # ind != 1
if ind != 1
n_para_rs += 1
end
end
Expand Down Expand Up @@ -224,7 +230,8 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr
initialize_cloud_settings!(cloud; tempered_update = tempered_update,
n_parts = n_parts, n_Φ = n_Φ, c = c, accept = target)

initialize_likelihoods!(loglikelihood, parameters, data, cloud; parallel = parallel)
initialize_likelihoods!(loglikelihood, parameters, data, cloud; parallel = parallel,
toggle = toggle)

elseif (tempered_update_prior_weight > 0.0) || (old_n_parts != n_parts)
# Resample from bridge distribution
Expand All @@ -242,7 +249,8 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr
# Make a cloud by drawing from the prior
if n_from_prior > 0
prior_cloud = Cloud(n_para, n_from_prior)
initial_draw!(loglikelihood, parameters, old_data, prior_cloud, parallel = parallel, regime_switching = regime_switching)
initial_draw!(loglikelihood, parameters, old_data, prior_cloud, parallel = parallel,
regime_switching = regime_switching, toggle = toggle)

cloud = Cloud(n_para, n_to_resample + n_from_prior)
update_cloud!(cloud, vcat(bridge_cloud.particles, prior_cloud.particles))
Expand All @@ -253,7 +261,8 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr
cloud = bridge_cloud
end

initialize_likelihoods!(loglikelihood, parameters, data, cloud; parallel = parallel)
initialize_likelihoods!(loglikelihood, parameters, data, cloud; parallel = parallel,
toggle = toggle)

reset_weights!(cloud)
cloud.resamples += 1
Expand All @@ -266,7 +275,8 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr
cloud = load(loadpath, "cloud")
else
# Instantiating Cloud object, update draws, loglh, & logpost
initial_draw!(loglikelihood, parameters, data, cloud; parallel = parallel, regime_switching = regime_switching)
initial_draw!(loglikelihood, parameters, data, cloud; parallel = parallel, regime_switching = regime_switching,
toggle = toggle)
initialize_cloud_settings!(cloud; tempered_update = tempered_update,
n_parts = n_parts, n_Φ = n_Φ, c = c, accept = target)
end
Expand Down Expand Up @@ -389,14 +399,16 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr
mutation_closure(cloud.particles[k, :], θ_bar_fr, R_fr, n_free_para,
blocks_free, blocks_all, ϕ_n, ϕ_n1; c = c, α = α,
n_mh_steps = n_mh_steps, old_data = old_data,
regime_switching = regime_switching)
regime_switching = regime_switching,
toggle = toggle)
end
else
hcat([mutation_closure(cloud.particles[k, :], θ_bar_fr, R_fr, n_free_para,
blocks_free, blocks_all, ϕ_n, ϕ_n1; c = c,
α = α, n_mh_steps = n_mh_steps,
old_data = old_data,
regime_switching = regime_switching) for k=1:n_parts]...)
regime_switching = regime_switching,
toggle = toggle) for k=1:n_parts]...)
end
update_cloud!(cloud, new_particles)
update_acceptance_rate!(cloud)
Expand Down
Loading

0 comments on commit b5ed5a4

Please sign in to comment.