-
Notifications
You must be signed in to change notification settings - Fork 10
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Error: SliceSampler ... shoud be initialized in the support #209
Comments
Thanks for the report and suggestions. It might be easier to give a recommendation after reading the model in question, potentially trying it locally as well. Conceptually, the slice sampler shrinks slides until they are within constraints, so we would expect that the once the sampler is in the support it should stay in it (hence the message talking about 'initialization' despite the problem occurring 'later' here). One possibility is that the prior sampler might propose configurations that have zero probability under the likelihood. Would that be possible in your case? If so, the best fix might be to restrict the prior to the set of parameters allowed by the likelihood function. PS: by default, the reason we did not do (a) and (b) by default is that they may void certain guarantees, e.g. correctness of log-normalization constant estimates. (b) might also lead to unbounded running time. However if we can identify a specific case where (a) would lead to correct result, I would be open to adding an option to facilitate it (potentially combined with warning messages). |
You can find the model here (just switch the NUTS sampler for pigeons and thats what i work with locally) if i understand correctly what you are saying, there might be an issue if the prior sampler makes a valid proposal but, in my case, the term in addlogprob is -Inf. so a fix (on my end) could be to catch infinite values and replace them with -1e300 and then hand them to addlogprob let me try that |
i guess the relevant part of my problem is this: Turing.@model function loglikelihood_function(data, m, observables)
z_ea ~ Truncated(InverseGamma(0.1, 2.0, μσ = true),0.01,3)
z_eb ~ Truncated(InverseGamma(0.1, 2.0, μσ = true),0.025,5)
...
kalman_prob = get_loglikelihood(m, data(observables), [z_ea, z_ea])
Turing.@addlogprob! kalman_prob
end your proposed fix doesn't work in my case because i don't know the parameter combinations which return a finite the fix I mentioned earlier ( ─────────────────────────────────
scans restarts Λ time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 2.27e-13 103 1.82e+11 1.14e-13 1 1 1 1
2 0 2.27e-13 103 1.82e+11 1.14e-13 1 1 1 1
ERROR: LoadError: ERROR: LoadError: Maximum number of iterations 1024 reached. Dumping info:
- Lbar = -5.1729024534740455
- Rbar = -5.172902453474045
- new_lp = -1.0e30
- z = -1.0e30
Stacktrace:
[1] error(s::String)
@ Base ./error.jl:35
[2] Maximum number of iterations 1024 reached. Dumping info:
- Lbar = -5.1729024534740455
- Rbar = -5.172902453474045
- new_lp = -1.0e30
- z = -1.0e30
Stacktrace:
[1] error(s::String)
@ Base ./error.jl:35
[2] slice_shrink!(h::P...) |
Thanks for the additional info, I'll have a look.. |
just some additional info: the above error came with 4 chains. with one chain i managed to let it run through |
Hi @thorek1 -- here it seems that you are trying to manually add a log-likelihood term in a Turing model. The proper way to do this is by making sure the Pigeons.jl/examples/turing-galaxy.jl Lines 21 to 33 in 764af10
Otherwise, the term would appear in both prior and posterior evaluations, making it impossible to apply tempering. |
so you suggest something like this!? Turing.@model function loglikelihood_function(data, m, observables)
z_ea ~ Truncated(InverseGamma(0.1, 2.0, μσ = true),0.01,3)
z_eb ~ Truncated(InverseGamma(0.1, 2.0, μσ = true),0.025,5)
...
if DynamicPPL.leafcontext(__context__) !== DynamicPPL.PriorContext()
kalman_prob = get_loglikelihood(m, data(observables), [z_ea, z_ea])
Turing.@addlogprob! kalman_prob
end
end let me give this a try |
Related to #206 |
now i remain stuck on a ───────────────────────────────────────────────────────────────────────────────────────
scans restarts time(s) allc(B) log(Z₁/Z₀) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 28.2 3.46e+10 0 NaN 1 1
2 0 28.2 3.46e+10 0 NaN 1 1
4 0 46.6 7.04e+10 0 NaN 1 1
4 0 46.6 7.04e+10 0 NaN 1 1
8 0 102 1.5e+11 0 NaN 1 1
8 0 102 1.5e+11 0 NaN 1 1
[164] signal (11.1): Segmentation fault
in expression starting at /home/cdsw/sample_sw07.jl:185
Allocations: 265281838 (Pool: 234497717; Big: 30784121); GC: 1940
[164] signal (11.1): Segmentation fault
in expression starting at /home/cdsw/sample_sw07.jl:185
Allocations: 265281838 (Pool: 234497717; Big: 30784121); GC: 1940
Segmentation fault (core dumped)
35584 any hint as to how to debug those? fix a seed and go step by step? |
Segfaults are the worst! Glad to hear there is progress though in terms of at least the preceding issue being resolved. A few questions/ideas:
To facilitate debugging, it might be helpful to use checkpoints; see https://pigeons.run/dev/checkpoints/ |
After fixing the context issue, I was able to run this model using multithreading (see my code at the end). Of course, due to computational limitations, I'm using just 4 chains and 4 rounds, but this shows that the model does work with Pigeons. In order to avoid the initialization in -Inf log-density region, I had to customize the initialization method. Basically, the idea is to sample from the prior until a non-infinite log-density state is found, and then use that as the initial point for all replicas. ============================================= # install latest versions of Pigeons and MacroModelling
using Pkg
Pkg.add([
PackageSpec(name="Pigeons", rev="main"),
PackageSpec(name="MacroModelling", rev="main")
])
using MacroModelling
import Turing
import Turing: NUTS, sample, logpdf, Truncated#, Normal, Beta, Gamma, InverseGamma,
using Random, CSV, DataFrames, Zygote, AxisKeys, MCMCChains
# using ComponentArrays, Optimization, OptimizationNLopt, OptimizationOptimisers
import DynamicPPL: DynamicPPL, logjoint
using Pigeons
@model SW07 begin
a[0] = calfa * rkf[0] + (1 - calfa) * (wf[0])
zcapf[0] = (1 / (czcap / (1 - czcap))) * rkf[0]
rkf[0] = wf[0] + labf[0] - kf[0]
kf[0] = kpf[-1] + zcapf[0]
invef[0] = (1 / (1 + cbetabar * cgamma)) * (invef[-1] + cbetabar * cgamma * invef[1] + (1 / (cgamma ^ 2 * csadjcost)) * pkf[0]) + qs[0]
pkf[0] = - rrf[0] + (1 / ((1 - chabb / cgamma) / (csigma * (1 + chabb / cgamma)))) * b[0] + (crk / (crk + (1 - ctou))) * rkf[1] + ((1 - ctou) / (crk + (1 - ctou))) * pkf[1]
cf[0] = (chabb / cgamma) / (1 + chabb / cgamma) * cf[-1] + (1 / (1 + chabb / cgamma)) * cf[1] + ((csigma - 1) * cwhlc / (csigma * (1 + chabb / cgamma))) * (labf[0] - labf[1]) - (1 - chabb / cgamma) / (csigma * (1 + chabb / cgamma)) * (rrf[0]) + b[0]
yf[0] = ccy * cf[0] + ciy * invef[0] + g[0] + crkky * zcapf[0]
yf[0] = cfc * (calfa * kf[0] + (1 - calfa) * labf[0] + a[0])
wf[0] = csigl * labf[0] + (1 / (1 - chabb / cgamma)) * cf[0] - (chabb / cgamma) / (1 - chabb / cgamma) * cf[-1]
kpf[0] = (1 - cikbar) * kpf[-1] + (cikbar) * invef[0] + (cikbar) * (cgamma ^ 2 * csadjcost) * qs[0]
mc[0] = calfa * rk[0] + (1 - calfa) * (w[0]) - a[0]
zcap[0] = (1 / (czcap / (1 - czcap))) * rk[0]
rk[0] = w[0] + lab[0] - k[0]
k[0] = kp[-1] + zcap[0]
inve[0] = (1 / (1 + cbetabar * cgamma)) * (inve[-1] + cbetabar * cgamma * inve[1] + (1 / (cgamma ^ 2 * csadjcost)) * pk[0]) + qs[0]
pk[0] = - r[0] + pinf[1] + (1 / ((1 - chabb / cgamma) / (csigma * (1 + chabb / cgamma)))) * b[0] + (crk / (crk + (1 - ctou))) * rk[1] + ((1 - ctou) / (crk + (1 - ctou))) * pk[1]
c[0] = (chabb / cgamma) / (1 + chabb / cgamma) * c[-1] + (1 / (1 + chabb / cgamma)) * c[1] + ((csigma - 1) * cwhlc / (csigma * (1 + chabb / cgamma))) * (lab[0] - lab[1]) - (1 - chabb / cgamma) / (csigma * (1 + chabb / cgamma)) * (r[0] - pinf[1]) + b[0]
y[0] = ccy * c[0] + ciy * inve[0] + g[0] + crkky * zcap[0]
y[0] = cfc * (calfa * k[0] + (1 - calfa) * lab[0] + a[0])
pinf[0] = (1 / (1 + cbetabar * cgamma * cindp)) * (cbetabar * cgamma * pinf[1] + cindp * pinf[-1] + ((1 - cprobp) * (1 - cbetabar * cgamma * cprobp) / cprobp) / ((cfc - 1) * curvp + 1) * (mc[0])) + spinf[0]
w[0] = (1 / (1 + cbetabar * cgamma)) * w[-1] + (cbetabar * cgamma / (1 + cbetabar * cgamma)) * w[1] + (cindw / (1 + cbetabar * cgamma)) * pinf[-1] - (1 + cbetabar * cgamma * cindw) / (1 + cbetabar * cgamma) * pinf[0] + (cbetabar * cgamma) / (1 + cbetabar * cgamma) * pinf[1] + (1 - cprobw) * (1 - cbetabar * cgamma * cprobw) / ((1 + cbetabar * cgamma) * cprobw) * (1 / ((clandaw - 1) * curvw + 1)) * (csigl * lab[0] + (1 / (1 - chabb / cgamma)) * c[0] - ((chabb / cgamma) / (1 - chabb / cgamma)) * c[-1] - w[0]) + sw[0]
r[0] = crpi * (1 - crr) * pinf[0] + cry * (1 - crr) * (y[0] - yf[0]) + crdy * (y[0] - yf[0] - y[-1] + yf[-1]) + crr * r[-1] + ms[0]
a[0] = crhoa * a[-1] + z_ea * ea[x]
b[0] = crhob * b[-1] + z_eb * eb[x]
g[0] = crhog * g[-1] + z_eg * eg[x] + cgy * z_ea * ea[x]
qs[0] = crhoqs * qs[-1] + z_eqs * eqs[x]
ms[0] = crhoms * ms[-1] + z_em * em[x]
spinf[0] = crhopinf * spinf[-1] + epinfma[0] - cmap * epinfma[-1]
epinfma[0] = z_epinf * epinf[x]
sw[0] = crhow * sw[-1] + ewma[0] - cmaw * ewma[-1]
ewma[0] = z_ew * ew[x]
kp[0] = (1 - cikbar) * kp[-1] + cikbar * inve[0] + cikbar * cgamma ^ 2 * csadjcost * qs[0]
dy[0] = y[0] - y[-1] + ctrend
dc[0] = c[0] - c[-1] + ctrend
dinve[0] = inve[0] - inve[-1] + ctrend
dw[0] = w[0] - w[-1] + ctrend
pinfobs[0] = (pinf[0]) + constepinf
robs[0] = (r[0]) + conster
labobs[0] = lab[0] + constelab
end
@parameters SW07 begin
ctou=.025
clandaw=1.5
cg=0.18
curvp=10
curvw=10
calfa=.24
# cgamma=1.004
# cbeta=.9995
csigma=1.5
# cpie=1.005
cfc=1.5
cgy=0.51
csadjcost= 6.0144
chabb= 0.6361
cprobw= 0.8087
csigl= 1.9423
cprobp= 0.6
cindw= 0.3243
cindp= 0.47
czcap= 0.2696
crpi= 1.488
crr= 0.8762
cry= 0.0593
crdy= 0.2347
crhoa= 0.9977
crhob= 0.5799
crhog= 0.9957
crhols= 0.9928
crhoqs= 0.7165
crhoas=1
crhoms=0
crhopinf=0
crhow=0
cmap = 0
cmaw = 0
clandap=cfc
cbetabar=cbeta*cgamma^(-csigma)
cr=cpie/(cbeta*cgamma^(-csigma))
crk=(cbeta^(-1))*(cgamma^csigma) - (1-ctou)
cw = (calfa^calfa*(1-calfa)^(1-calfa)/(clandap*crk^calfa))^(1/(1-calfa))
cikbar=(1-(1-ctou)/cgamma)
cik=(1-(1-ctou)/cgamma)*cgamma
clk=((1-calfa)/calfa)*(crk/cw)
cky=cfc*(clk)^(calfa-1)
ciy=cik*cky
ccy=1-cg-cik*cky
crkky=crk*cky
cwhlc=(1/clandaw)*(1-calfa)/calfa*crk*cky/ccy
cwly=1-crk*cky
conster=(cr-1)*100
# ctrend=(cgamma-1)*100
ctrend=(1.004-1)*100
# constepinf=(cpie-1)*100
constepinf=(1.005-1)*100
cpie=1+constepinf/100
cgamma=1+ctrend/100
cbeta=1/(1+constebeta/100)
constebeta = 100 / .9995 - 100
constelab=0
z_ea = 0.4618
z_eb = 1.8513
z_eg = 0.6090
z_eqs = 0.6017
z_em = 0.2397
z_epinf = 0.1455
z_ew = 0.2089
end
# load data
dat = CSV.read(joinpath(pkgdir(MacroModelling), "test", "data", "usmodel.csv"), DataFrame)
data = KeyedArray(Array(dat)',Variable = Symbol.(strip.(names(dat))), Time = 1:size(dat)[1])
# declare observables
observables = [:dy, :dc, :dinve, :labobs, :pinfobs, :dw, :robs]
# Subsample from 1966Q1 - 2004Q4
# subset observables in data
data = data(observables,75:230)
# functions to map mean and standard deviations to distribution parameters
DynamicPPL.@model function SW07_loglikelihood_function(data, m, observables,fixed_parameters)
z_ea ~ Truncated(InverseGamma(0.1, 2.0, μσ = true),0.01,3)
z_eb ~ Truncated(InverseGamma(0.1, 2.0, μσ = true),0.025,5)
z_eg ~ Truncated(InverseGamma(0.1, 2.0, μσ = true),0.01,3)
z_eqs ~ Truncated(InverseGamma(0.1, 2.0, μσ = true),0.01,3)
z_em ~ Truncated(InverseGamma(0.1, 2.0, μσ = true),0.01,3)
z_epinf ~ Truncated(InverseGamma(0.1, 2.0, μσ = true),0.01,3)
z_ew ~ Truncated(InverseGamma(0.1, 2.0, μσ = true),0.01,3)
crhoa ~ Truncated(Beta(0.5, 0.20, μσ = true),.01,.9999)
crhob ~ Truncated(Beta(0.5, 0.20, μσ = true),.01,.9999)
crhog ~ Truncated(Beta(0.5, 0.20, μσ = true),.01,.9999)
crhoqs ~ Truncated(Beta(0.5, 0.20, μσ = true),.01,.9999)
crhoms ~ Truncated(Beta(0.5, 0.20, μσ = true),.01,.9999)
crhopinf~ Truncated(Beta(0.5, 0.20, μσ = true),.01,.9999)
crhow ~ Truncated(Beta(0.5, 0.20, μσ = true),.001,.9999)
cmap ~ Truncated(Beta(0.5, 0.2, μσ = true),0.01,.9999)
cmaw ~ Truncated(Beta(0.5, 0.2, μσ = true),0.01,.9999)
csadjcost~ Truncated(Normal(4,1.5),2,15)
csigma ~ Truncated(Normal(1.50,0.375),0.25,3)
chabb ~ Truncated(Beta(0.7, 0.1, μσ = true),0.001,0.99)
cprobw ~ Truncated(Beta(0.5, 0.1, μσ = true),0.3,0.95)
csigl ~ Truncated(Normal(2,0.75),0.25,10)
cprobp ~ Truncated(Beta(0.5, 0.10, μσ = true),0.5,0.95)
cindw ~ Truncated(Beta(0.5, 0.15, μσ = true),0.01,0.99)
cindp ~ Truncated(Beta(0.5, 0.15, μσ = true),0.01,0.99)
czcap ~ Truncated(Beta(0.5, 0.15, μσ = true),0.01,1)
cfc ~ Truncated(Normal(1.25,0.125),1.0,3)
crpi ~ Truncated(Normal(1.5,0.25),1.0,3)
crr ~ Truncated(Beta(0.75, 0.10, μσ = true),0.5,0.975)
cry ~ Truncated(Normal(0.125,0.05),0.001,0.5)
crdy ~ Truncated(Normal(0.125,0.05),0.001,0.5)
constepinf~ Truncated(Gamma(0.625,0.1, μσ = true),0.1,2.0)
constebeta~ Truncated(Gamma(0.25,0.1, μσ = true),0.01,2.0)
constelab ~ Truncated(Normal(0.0,2.0),-10.0,10.0)
ctrend ~ Truncated(Normal(0.4,0.10),0.1,0.8)
cgy ~ Truncated(Normal(0.5,0.25),0.01,2.0)
calfa ~ Truncated(Normal(0.3,0.05),0.01,1.0)
ctou, clandaw, cg, curvp, curvw, crhols, crhoas = fixed_parameters
parameters_combined = [ctou,clandaw,cg,curvp,curvw,calfa,csigma,cfc,cgy,csadjcost,chabb,cprobw,csigl,cprobp,cindw,cindp,czcap,crpi,crr,cry,crdy,crhoa,crhob,crhog,crhols,crhoqs,crhoas,crhoms,crhopinf,crhow,cmap,cmaw,constelab,z_ea,z_eb,z_eg,z_eqs,z_em,z_epinf,z_ew,ctrend,constepinf,constebeta]
if DynamicPPL.leafcontext(__context__) !== DynamicPPL.PriorContext()
kalman_prob = get_loglikelihood(m, data(observables), parameters_combined)
# println(kalman_prob)
DynamicPPL.@addlogprob! kalman_prob
end
end
SW07.parameter_values[indexin([:crhoms, :crhopinf, :crhow, :cmap, :cmaw],SW07.parameters)] .= 0.02
fixed_parameters = SW07.parameter_values[indexin([:ctou,:clandaw,:cg,:curvp,:curvw,:crhols,:crhoas],SW07.parameters)]
SW07_loglikelihood = SW07_loglikelihood_function(data, SW07, observables, fixed_parameters)
###############################################################################
# inference using Pigeons
###############################################################################
# generate a Pigeons log potential
sw07_lp = TuringLogPotential(SW07_loglikelihood)
# find a feasible starting point
pt = pigeons(target = sw07_lp, n_rounds = 0, n_chains = 1)
replica = pt.replicas[end]
xmax = deepcopy(replica.state)
lpmax = sw07_lp(xmax)
i = 0
while !isfinite(lpmax) && i < 100
Pigeons.sample_iid!(sw07_lp, replica, pt.shared)
new_lp = sw07_lp(replica.state)
if new_lp > lpmax
lpmax = new_lp
xmax = deepcopy(replica.state)
end
i += 1
end
# define a specific initialization for this model
Pigeons.initialization(::TuringLogPotential{typeof(SW07_loglikelihood)}, ::AbstractRNG, ::Int64) = deepcopy(xmax)
# run Pigeons
pt = pigeons(target = sw07_lp, n_chains = 4, n_rounds = 4, multithreaded=true) |
neat. thanks a lot! i will let it run on a couple of use cases over night. any chance of not having to go that route in the future? as in, not using the context code within the turing model and not needing to find the feasible point first? |
No worries, glad I could help. Regarding your questions:
|
Nice!! Creating a new issue for 2... (since it is a bit more general than slice sampling / Turing) |
so far the proposed fixes work, rarely, as in most of the times i get segfaults. I am using checkpoints, am setting a seed, and will, time allowing, investigate why this still fails (on linux, with chains>1 and multithreading). any hints or best practices are most welcome |
Sorry to hear. Feel free to send us instructions to see if we can replicate these segfaults for investigation. Based on an earlier message, I was under the impression that segfault was occurring even in the single thread context. Is this still true or they only occur in a multi-thread context? |
collecting evidence... checking without multithreading now |
here is what i am running (with this data file: usmodel.csv): using MKL
# using LinearAlgebra
# BLAS.get_config()
using MacroModelling
using Serialization
using StatsPlots
import Turing
import Turing: NUTS, PG, IS, sample, logpdf, Truncated#, Normal, Beta, Gamma, InverseGamma,
using Random, CSV, DataFrames, Zygote, AxisKeys, MCMCChains
# using ComponentArrays, Optimization, OptimizationNLopt, OptimizationOptimisers
import DynamicPPL: logjoint
import DynamicPPL
import ChainRulesCore: @ignore_derivatives, ignore_derivatives
import Pigeons
import Random
Random.seed!(1)
println("Threads used: ", Threads.nthreads())
smpler = "pigeons" #
smple = "original" #
mdl = "linear" #
chns = 1 #
if smple == "original"
smpl = "1966Q1-2004Q4"
sample_idx = 75:230
dat = CSV.read("usmodel.csv", DataFrame)
end
# define callback
# Define the path for the CSV file
csv_file_path = "sw07_$(mdl)_$(smpler)_$(smpl)_samples.csv"
# Initialize a DataFrame to store the data
df = DataFrame(iteration = Float64[])
function callback(rng, model, sampler, sample, state, i; kwargs...)
# Prepare a row for the DataFrame
row = Dict("iteration" => Float64(i))
for (name, value) in sample.θ
row[string(name)] = value
end
# If the DataFrame `df` does not have columns for these names, add them
for name in keys(row)
if !any(name .== names(df))
df[!, name] = Union{Missing, Any}[missing for _ in 1:nrow(df)]
end
end
# Append the new data to the DataFrame
push!(df, row)
# Write the updated DataFrame to the CSV file
# Note: To avoid performance issues, consider writing periodically instead of on every callback
CSV.write(csv_file_path, df, append=true)
end
@model SW07 begin
a[0] = calfa * rkf[0] + (1 - calfa) * (wf[0])
zcapf[0] = (1 / (czcap / (1 - czcap))) * rkf[0]
rkf[0] = wf[0] + labf[0] - kf[0]
kf[0] = kpf[-1] + zcapf[0]
invef[0] = (1 / (1 + cbetabar * cgamma)) * (invef[-1] + cbetabar * cgamma * invef[1] + (1 / (cgamma ^ 2 * csadjcost)) * pkf[0]) + qs[0]
pkf[0] = - rrf[0] + (1 / ((1 - chabb / cgamma) / (csigma * (1 + chabb / cgamma)))) * b[0] + (crk / (crk + (1 - ctou))) * rkf[1] + ((1 - ctou) / (crk + (1 - ctou))) * pkf[1]
cf[0] = (chabb / cgamma) / (1 + chabb / cgamma) * cf[-1] + (1 / (1 + chabb / cgamma)) * cf[1] + ((csigma - 1) * cwhlc / (csigma * (1 + chabb / cgamma))) * (labf[0] - labf[1]) - (1 - chabb / cgamma) / (csigma * (1 + chabb / cgamma)) * (rrf[0]) + b[0]
yf[0] = ccy * cf[0] + ciy * invef[0] + g[0] + crkky * zcapf[0]
yf[0] = cfc * (calfa * kf[0] + (1 - calfa) * labf[0] + a[0])
wf[0] = csigl * labf[0] + (1 / (1 - chabb / cgamma)) * cf[0] - (chabb / cgamma) / (1 - chabb / cgamma) * cf[-1]
kpf[0] = (1 - cikbar) * kpf[-1] + (cikbar) * invef[0] + (cikbar) * (cgamma ^ 2 * csadjcost) * qs[0]
mc[0] = calfa * rk[0] + (1 - calfa) * (w[0]) - a[0]
zcap[0] = (1 / (czcap / (1 - czcap))) * rk[0]
rk[0] = w[0] + lab[0] - k[0]
k[0] = kp[-1] + zcap[0]
inve[0] = (1 / (1 + cbetabar * cgamma)) * (inve[-1] + cbetabar * cgamma * inve[1] + (1 / (cgamma ^ 2 * csadjcost)) * pk[0]) + qs[0]
pk[0] = - r[0] + pinf[1] + (1 / ((1 - chabb / cgamma) / (csigma * (1 + chabb / cgamma)))) * b[0] + (crk / (crk + (1 - ctou))) * rk[1] + ((1 - ctou) / (crk + (1 - ctou))) * pk[1]
c[0] = (chabb / cgamma) / (1 + chabb / cgamma) * c[-1] + (1 / (1 + chabb / cgamma)) * c[1] + ((csigma - 1) * cwhlc / (csigma * (1 + chabb / cgamma))) * (lab[0] - lab[1]) - (1 - chabb / cgamma) / (csigma * (1 + chabb / cgamma)) * (r[0] - pinf[1]) + b[0]
y[0] = ccy * c[0] + ciy * inve[0] + g[0] + crkky * zcap[0]
y[0] = cfc * (calfa * k[0] + (1 - calfa) * lab[0] + a[0])
pinf[0] = (1 / (1 + cbetabar * cgamma * cindp)) * (cbetabar * cgamma * pinf[1] + cindp * pinf[-1] + ((1 - cprobp) * (1 - cbetabar * cgamma * cprobp) / cprobp) / ((cfc - 1) * curvp + 1) * (mc[0])) + spinf[0]
w[0] = (1 / (1 + cbetabar * cgamma)) * w[-1] + (cbetabar * cgamma / (1 + cbetabar * cgamma)) * w[1] + (cindw / (1 + cbetabar * cgamma)) * pinf[-1] - (1 + cbetabar * cgamma * cindw) / (1 + cbetabar * cgamma) * pinf[0] + (cbetabar * cgamma) / (1 + cbetabar * cgamma) * pinf[1] + (1 - cprobw) * (1 - cbetabar * cgamma * cprobw) / ((1 + cbetabar * cgamma) * cprobw) * (1 / ((clandaw - 1) * curvw + 1)) * (csigl * lab[0] + (1 / (1 - chabb / cgamma)) * c[0] - ((chabb / cgamma) / (1 - chabb / cgamma)) * c[-1] - w[0]) + sw[0]
r[0] = crpi * (1 - crr) * pinf[0] + cry * (1 - crr) * (y[0] - yf[0]) + crdy * (y[0] - yf[0] - y[-1] + yf[-1]) + crr * r[-1] + ms[0]
a[0] = crhoa * a[-1] + z_ea * ea[x]
b[0] = crhob * b[-1] + z_eb * eb[x]
g[0] = crhog * g[-1] + z_eg * eg[x] + cgy * z_ea * ea[x]
qs[0] = crhoqs * qs[-1] + z_eqs * eqs[x]
ms[0] = crhoms * ms[-1] + z_em * em[x]
spinf[0] = crhopinf * spinf[-1] + epinfma[0] - cmap * epinfma[-1]
epinfma[0] = z_epinf * epinf[x]
sw[0] = crhow * sw[-1] + ewma[0] - cmaw * ewma[-1]
ewma[0] = z_ew * ew[x]
kp[0] = (1 - cikbar) * kp[-1] + cikbar * inve[0] + cikbar * cgamma ^ 2 * csadjcost * qs[0]
dy[0] = y[0] - y[-1] + ctrend
dc[0] = c[0] - c[-1] + ctrend
dinve[0] = inve[0] - inve[-1] + ctrend
dw[0] = w[0] - w[-1] + ctrend
pinfobs[0] = (pinf[0]) + constepinf
robs[0] = (r[0]) + conster
labobs[0] = lab[0] + constelab
end
@parameters SW07 begin
ctou=.025
clandaw=1.5
cg=0.18
curvp=10
curvw=10
calfa=.24
# cgamma=1.004
# cbeta=.9995
csigma=1.5
# cpie=1.005
cfc=1.5
cgy=0.51
csadjcost= 6.0144
chabb= 0.6361
cprobw= 0.8087
csigl= 1.9423
cprobp= 0.6
cindw= 0.3243
cindp= 0.47
czcap= 0.2696
crpi= 1.488
crr= 0.8762
cry= 0.0593
crdy= 0.2347
crhoa= 0.9977
crhob= 0.5799
crhog= 0.9957
crhols= 0.9928
crhoqs= 0.7165
crhoas=1
crhoms=0
crhopinf=0
crhow=0
cmap = 0
cmaw = 0
clandap=cfc
cbetabar=cbeta*cgamma^(-csigma)
cr=cpie/(cbeta*cgamma^(-csigma))
crk=(cbeta^(-1))*(cgamma^csigma) - (1-ctou)
cw = (calfa^calfa*(1-calfa)^(1-calfa)/(clandap*crk^calfa))^(1/(1-calfa))
cikbar=(1-(1-ctou)/cgamma)
cik=(1-(1-ctou)/cgamma)*cgamma
clk=((1-calfa)/calfa)*(crk/cw)
cky=cfc*(clk)^(calfa-1)
ciy=cik*cky
ccy=1-cg-cik*cky
crkky=crk*cky
cwhlc=(1/clandaw)*(1-calfa)/calfa*crk*cky/ccy
cwly=1-crk*cky
conster=(cr-1)*100
# ctrend=(cgamma-1)*100
ctrend=(1.004-1)*100
# constepinf=(cpie-1)*100
constepinf=(1.005-1)*100
cpie=1+constepinf/100
cgamma=1+ctrend/100
cbeta=1/(1+constebeta/100)
constebeta = 100 / .9995 - 100
constelab=0
z_ea = 0.4618
z_eb = 1.8513
z_eg = 0.6090
z_eqs = 0.6017
z_em = 0.2397
z_epinf = 0.1455
z_ew = 0.2089
end
# load data
data = KeyedArray(Array(dat)',Variable = Symbol.(strip.(names(dat))), Time = 1:size(dat)[1])
# declare observables
observables = [:dy, :dc, :dinve, :labobs, :pinfobs, :dw, :robs]
# Subsample
# subset observables in data
data = data(observables, sample_idx)
# functions to map mean and standard deviations to distribution parameters
Turing.@model function SW07_loglikelihood_function(data, m, observables,fixed_parameters)
z_ea ~ Truncated(InverseGamma(0.1, 2.0, μσ = true),0.01,3)
z_eb ~ Truncated(InverseGamma(0.1, 2.0, μσ = true),0.025,5)
z_eg ~ Truncated(InverseGamma(0.1, 2.0, μσ = true),0.01,3)
z_eqs ~ Truncated(InverseGamma(0.1, 2.0, μσ = true),0.01,3)
z_em ~ Truncated(InverseGamma(0.1, 2.0, μσ = true),0.01,3)
z_epinf ~ Truncated(InverseGamma(0.1, 2.0, μσ = true),0.01,3)
z_ew ~ Truncated(InverseGamma(0.1, 2.0, μσ = true),0.01,3)
crhoa ~ Truncated(Beta(0.5, 0.20, μσ = true),.01,.9999)
crhob ~ Truncated(Beta(0.5, 0.20, μσ = true),.01,.9999)
crhog ~ Truncated(Beta(0.5, 0.20, μσ = true),.01,.9999)
crhoqs ~ Truncated(Beta(0.5, 0.20, μσ = true),.01,.9999)
crhoms ~ Truncated(Beta(0.5, 0.20, μσ = true),.01,.9999)
crhopinf~ Truncated(Beta(0.5, 0.20, μσ = true),.01,.9999)
crhow ~ Truncated(Beta(0.5, 0.20, μσ = true),.001,.9999)
cmap ~ Truncated(Beta(0.5, 0.2, μσ = true),0.01,.9999)
cmaw ~ Truncated(Beta(0.5, 0.2, μσ = true),0.01,.9999)
csadjcost~ Truncated(Normal(4,1.5),2,15)
csigma ~ Truncated(Normal(1.50,0.375),0.25,3)
chabb ~ Truncated(Beta(0.7, 0.1, μσ = true),0.001,0.99)
cprobw ~ Truncated(Beta(0.5, 0.1, μσ = true),0.3,0.95)
csigl ~ Truncated(Normal(2,0.75),0.25,10)
cprobp ~ Truncated(Beta(0.5, 0.10, μσ = true),0.5,0.95)
cindw ~ Truncated(Beta(0.5, 0.15, μσ = true),0.01,0.99)
cindp ~ Truncated(Beta(0.5, 0.15, μσ = true),0.01,0.99)
czcap ~ Truncated(Beta(0.5, 0.15, μσ = true),0.01,1)
cfc ~ Truncated(Normal(1.25,0.125),1.0,3)
crpi ~ Truncated(Normal(1.5,0.25),1.0,3)
crr ~ Truncated(Beta(0.75, 0.10, μσ = true),0.5,0.975)
cry ~ Truncated(Normal(0.125,0.05),0.001,0.5)
crdy ~ Truncated(Normal(0.125,0.05),0.001,0.5)
constepinf~ Truncated(Gamma(0.625,0.1, μσ = true),0.1,2.0)
constebeta~ Truncated(Gamma(0.25,0.1, μσ = true),0.01,2.0)
constelab ~ Truncated(Normal(0.0,2.0),-10.0,10.0)
ctrend ~ Truncated(Normal(0.4,0.10),0.1,0.8)
cgy ~ Truncated(Normal(0.5,0.25),0.01,2.0)
calfa ~ Truncated(Normal(0.3,0.05),0.01,1.0)
if DynamicPPL.leafcontext(__context__) !== DynamicPPL.PriorContext()
ctou, clandaw, cg, curvp, curvw, crhols, crhoas = fixed_parameters
parameters_combined = [ctou,clandaw,cg,curvp,curvw,calfa,csigma,cfc,cgy,csadjcost,chabb,cprobw,csigl,cprobp,cindw,cindp,czcap,crpi,crr,cry,crdy,crhoa,crhob,crhog,crhols,crhoqs,crhoas,crhoms,crhopinf,crhow,cmap,cmaw,constelab,z_ea,z_eb,z_eg,z_eqs,z_em,z_epinf,z_ew,ctrend,constepinf,constebeta]
kalman_prob = get_loglikelihood(m, data(observables), parameters_combined)
Turing.@addlogprob! kalman_prob
end
end
SW07.parameter_values[indexin([:crhoms, :crhopinf, :crhow, :cmap, :cmaw],SW07.parameters)] .= 0.02
fixed_parameters = SW07.parameter_values[indexin([:ctou,:clandaw,:cg,:curvp,:curvw,:crhols,:crhoas],SW07.parameters)]
dir_name = "sw07_$(mdl)_$(smpler)_$(smpl)_samples_$(chns)_chains"
if !isdir(dir_name) mkdir(dir_name) end
cd(dir_name)
println("Current working directory: ", pwd())
SW07_loglikelihood = SW07_loglikelihood_function(data, SW07, observables, fixed_parameters)
if smpler == "pigeons"
# generate a Pigeons log potential
sw07_lp = Pigeons.TuringLogPotential(SW07_loglikelihood)
# find a feasible starting point
pt = Pigeons.pigeons(target = sw07_lp, n_rounds = 0, n_chains = 1)
replica = pt.replicas[end]
XMAX = deepcopy(replica.state)
LPmax = sw07_lp(XMAX)
i = 0
while !isfinite(LPmax) && i < 1000
Pigeons.sample_iid!(sw07_lp, replica, pt.shared)
new_LP = sw07_lp(replica.state)
if new_LP > LPmax
LPmax = new_LP
XMAX = deepcopy(replica.state)
end
i += 1
end
# define a specific initialization for this model
Pigeons.initialization(::Pigeons.TuringLogPotential{typeof(SW07_loglikelihood)}, ::AbstractRNG, ::Int64) = deepcopy(XMAX)
pt = Pigeons.pigeons(target = sw07_lp,
checkpoint = true,
record = [Pigeons.traces; Pigeons.round_trip; Pigeons.record_default(); Pigeons.disk],
n_chains = chns,
n_rounds = 10,
multithreaded = false)
samps = MCMCChains.Chains(Pigeons.get_sample(pt))
end
serialize("samples.jls", samps)
my_plot = StatsPlots.plot(samps)
StatsPlots.savefig(my_plot, "samples.png") |
Thank you! The manifest.toml file would be great too, to make sure we use same versions. |
there you go
|
and i can now confirm single threaded runs also fail with segfaults and are not reproducible with the seed (they fail at different rounds) |
I am using a Turing model and it can happen especially at later stages/scans or when i use >1 chains that I get the below error (deleted parts to make it fit here).
how can I avoid this?
it appears he proposes a sample outside the prior bounds which returns infinite log probability -> error.
I dont know what the best solution is but from a user perspective it would be great if sampling continues. as a proposal: would it be possible to either a) ignore the error, or b) take draws until you get finite log probability?
The text was updated successfully, but these errors were encountered: