Lbar becoming larger than Rbar in SliceSampler causing max iterations reached #245
Comments
|
Hi -- that output means that the slicer is stuck trying to find a point with log probability higher than 10^26, which is basically inf. In turn, this suggests that there may be an issue with the log potential. Increasing the iterations likely won't help. It is better to understand why there exist a point with seemingly infinite log potential. If you post a reproducible example we might be able to help. |
|
@miguelbiron can you make a PR to improve that error message? It is certainly confusing if it only says "max iterations reached" -- not clear what "iterations" it's talking about here. Probably want to say something like "slicer maxed out iters while trying to ....". Also want to include a helpful hint about what is going on and how one might fix. |
|
@trevorcampbell ok that sounds good! |
|
Looking at the Blang implementation I see I had observed cases where due to precision errors the shrinkage can fail, and so I had a check that if the left and right end point are approximately equal, just stay at the old point. (https://github.com/UBC-Stat-ML/blangSDK/blob/b8642c9c2a0adab8a5b6da96f2a7889f1b81b6cc/src/main/java/blang/mcmc/RealSliceSampler.java#L111) I'll do a PR to do the same thing in Pigeons unless there are objections... |
alexandrebouchard
changed the title
Sorry I missed this. using CSV
using DataFrames
using Turing
using Pigeons
using SequentialSamplingModels
using Downloads
using LinearAlgebra
using Statistics
function data_poly(x, degree=2; orthogonal=false)
if orthogonal
z = x .- mean(x) # Center the data by subtracting its mean
X = hcat([z .^ deg for deg in 1:degree]...) # Create the matrix of powers up to 'degree'
QR = qr(X) # Perform QR decomposition
X = Matrix(QR.Q) # Extract the orthogonal matrix Q
else
X = hcat([x .^ deg for deg in 1:degree]...) # Create the matrix of powers up to 'degree'
end
return X
end
@model function model_exgaussian(data; min_rt=minimum(data.rt), isi=nothing)
# Transform ISI into polynomials
isi = data_poly(isi, 2; orthogonal=true)
# Priors for coefficients
drift_intercept ~ Normal(0, 1)
drift_isi1 ~ Normal(0, 1)
drift_isi2 ~ Normal(0, 1)
σ ~ Normal(0, 1)
τ ~ Normal(log(0.2), 1)
for i in 1:length(data)
drift = drift_intercept
drift += drift_isi1 * isi[i, 1]
drift += drift_isi2 * isi[i, 2]
data[i] ~ ExGaussian(exp(drift), exp(σ), exp(τ))
end
end
df = CSV.read(Downloads.download("https://raw.githubusercontent.com/RealityBending/DoggoNogo/main/study1/data/data_game.csv"), DataFrame)
fit = model_exgaussian(df.RT, min_rt=minimum(df.RT), isi=df.ISI)
pt = pigeons(target=TuringLogPotential(fit);
record=[Pigeons.traces],
n_rounds=5,
n_chains=10,
# checkpoint=true,
multithreaded=true,
seed=123)I'm fairly convinced that the model is tightly parametrized. @itsdfish could you correct me if I'm wrong? |
|
Good news, the PR I created to address this issue, PR #268, appears to fix that issue. Running your example I now get: Thank you so much for pointing this out! Relatedly, we have started working on a collection of difficult posterior distributions to help test and improve pigeons. (see beta version here). Would you be interested in having the above model part of this collection? If so, we would just need a bib reference entry to provide credits. |
|
Using the latest However, it often spat this new message And the traces look funny:
MWE: using CSV
using DataFrames
using Turing
using Pigeons
using SequentialSamplingModels
using Downloads
using LinearAlgebra
using Statistics
using StatsPlots
function data_poly(x, degree=2; orthogonal=false)
if orthogonal
z = x .- mean(x) # Center the data by subtracting its mean
X = hcat([z .^ deg for deg in 1:degree]...) # Create the matrix of powers up to 'degree'
QR = qr(X) # Perform QR decomposition
X = Matrix(QR.Q) # Extract the orthogonal matrix Q
else
X = hcat([x .^ deg for deg in 1:degree]...) # Create the matrix of powers up to 'degree'
end
return X
end
@model function model_ExGaussian(rt; min_rt=minimum(data.rt), isi=nothing)
# Transform ISI into polynomials
isi = data_poly(isi, 2; orthogonal=true)
# Priors for coefficients
μ_intercept ~ Normal(0, 1)
μ_isi1 ~ Normal(0, 1)
μ_isi2 ~ Normal(0, 1)
σ ~ Normal(0, 1)
τ ~ Normal(log(0.2), 1)
for i in 1:length(rt)
μ = μ_intercept
μ += μ_isi1 * isi[i, 1]
μ += μ_isi2 * isi[i, 2]
rt[i] ~ ExGaussian(μ, exp(σ), exp(τ))
end
end
df = CSV.read(Downloads.download("https://raw.githubusercontent.com/RealityBending/DoggoNogo/main/study1/data/data_game.csv"), DataFrame)
fit = model_ExGaussian(df.RT[1:100, :], min_rt=minimum(df.RT), isi=df.ISI[1:100, :])
pt = pigeons(target=TuringLogPotential(fit);
record=[Pigeons.traces],
n_rounds=9,
n_chains=10,
multithreaded=true,
seed=123)
posteriors = Chains(pt)
StatsPlots.plot(posteriors; size=(800, 800))Any pointers as to what I'm doing wrong? |
|
For the "set of report has changed", this should be fixed in commit 1e0484f which I pushed earlier today. The underlying phenomenon is the following:
It's helpful to know which of these 2 you are facing. Such target is intrinsically difficult. Contrary to popular belief, HMC often struggle on those as well. |
|
To find which of the 2 you are facing, you can try running on a small subset of the data... |
|
@miguelbiron, good question. Here is the logpdf function. I'm not quite sure if anything is problematic. The function returns -Inf when tau=0 because of division by zero and the exponential component of the exgaussian is undefined at zero. Does anything look like an issue to you? One potential issue might be numberical instability here |
|
I truncated the priors for sigma and tau: I had success with NUTS: The ratio log(z_1 / z_0) was very large for Pigeons.jl. I suspect this is indicative of a problem. |
|
Apologies for the bump, but is it possible that #273 will also help this case or is it completely unrelated? |
|
Hi @DominiqueMakowski -- I don't know if @itsdfish has ruled out yet numerical instabilities as the issues here. But if they are not -- and I have a feeling that they might not be at least the main culprit, given my comment here -- then it is possible that the models from this package tend to have likelihoods with stricter smaller support than the prior. If that's the case, then you don't need to wait for #273 to test that; you can just define a custom deterministic initialization for your particular target like so: fit = model_ExGaussian(df.RT[1:100, :], min_rt=minimum(df.RT), isi=df.ISI[1:100, :])
target = TuringLogPotential(fit)
Pigeons.initialization(::typeof(target), ::AbstractRNG, ::Int64) = zeros(DIMENSION)where you would need to manually provide the dimension. It does not have to be just zeros either. Just remember that these values are assumed to be in unconstrained space. |
|
Oh but you know what, your issue is occuring at a different point, after the initialization check so disregard the previous comment. |
|
You can also try switching the explorer to AutoMALA. |
|
I can confirm that adding (note that using Enzyme as in the docs fails for my usecase but hopefully this will get fixed as Enzyme continues to be integrated) |
|
I'm glad it worked! Yes, enzyme is not working yet in Turing, that's a work in progress on their end. You may get better performance with ReverseDiff, but make sure that you get the same results as woth forwarddiff first. |
|
Damn it I rejoiced too fast, it warked on a subset of the data but on the whole data it fails with: I'll try the initialization trick, but could you clarify what you mean by : where you would need to manually provide the dimension what does the dimension correspond to? |
|
Its the total number of parameters that you're inferring. |
|
Ah no wait the code I gave above won't work as is, it requires a special construction for turing models. I'm busy now but I can provide an example during the day. |
|
Ok this should work fit = model_ExGaussian(df.RT[1:100, :], min_rt=minimum(df.RT), isi=df.ISI[1:100, :])
fit_target = TuringLogPotential(fit)
function Pigeons.initialization(target::typeof(fit_target), rng::AbstractRNG, ::Int64)
vi = DynamicPPL.VarInfo(rng, target.model, DynamicPPL.SampleFromPrior(), DynamicPPL.PriorContext())
DynamicPPL.link!!(vi, DynamicPPL.SampleFromPrior(), target.model)
d = length(DynamicPPL.getall(vi)) # get dimension
new_vals = zeros(d)
DynamicPPL.setall!(vi, new_vals)
lp = DynamicPPL.logjoint(target.model, vi) # recompute logjoint with new values
@assert isfinite(lp)
return vi
end |
|
Unfortunately despite adding that function it fails with a new error after some time: Note that the same model gets sampled nicely within a few seconds using regular MCMC so again it does seem to me like it's a model specification issue... |
I've been having issues with running
pigeons()on a Turing model and the following error keeps popping up:How do we increase the number of iterations?
I tried adding
max_iterargument to thepigeons()call but it seems like it doesn't workthanks!
The text was updated successfully, but these errors were encountered: