Unexpected AssertionError

[astrozot]

I apologize, this is more a help request than a bug report or an enhancement suggestion.

During a Bayesian inference of parameters of a fairly complex log-posterior I encountered a problem that I do not really understand:

─────────────────────────────────────────────────────────────────────────────────────────────────────────────
  scans     restarts      Λ        time(s)    allc(B)  log(Z₁/Z₀)   min(α)     mean(α)    min(αₑ)   mean(αₑ) 
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
        2          0          4      0.224   7.08e+07  -4.86e+04          0      0.556          1          1 
        4          0          4      0.808   2.55e+08  -2.33e+04          0      0.556          1          1 
        8          0        7.8       1.71   5.97e+08  -1.66e+04          0      0.133          1          1 
       16          0       8.97       3.97   1.25e+09  -1.33e+04          0    0.00376          1          1 
       32          0       8.86       8.99   2.76e+09  -1.26e+04          0     0.0156          1          1 
ERROR: AssertionError: isfinite(z)
Stacktrace:
  [1] slice_shrink!(h::Pigeons.SliceSampler, replica::Pigeons.Replica{Vector{Float64}, NamedTuple{(:explorer_acceptance_pr, :explorer_n_steps, :swap_acceptance_pr, :log_sum_ratio, :traces, :index_process, :round_trip, :timing_extrema, :allocation_extrema, :_transformed_online), Tuple{OnlineStatsBase.GroupBy{Int64, Number, OnlineStatsBase.Mean{Float64, OnlineStatsBase.EqualWeight}}, OnlineStatsBase.GroupBy{Int64, Number, OnlineStatsBase.Sum{Float64}}, OnlineStatsBase.GroupBy{Tuple{Int64, Int64}, Number, OnlineStatsBase.Mean{Float64, OnlineStatsBase.EqualWeight}}, OnlineStatsBase.GroupBy{Tuple{Int64, Int64}, Number, Pigeons.LogSum{Float64}}, Dict{Pair{Int64, Int64}, Any}, Dict{Int64, Vector{Int64}}, Pigeons.RoundTripRecorder, NonReproducible{...}, NonReproducible{...}, Pigeons.OnlineStateRecorder}}}, z::Float64, L::Float64, R::Float64, lp_L::Float64, lp_R::Float64, pointer::Base.RefArray{Float64, Vector{Float64}, Nothing}, log_potential::InterpolatedLogPotential{...})
    @ Pigeons ~/.julia/packages/Pigeons/rBo9q/src/explorers/SliceSampler.jl:153
  [2] slice_sample_coord!(h::Pigeons.SliceSampler, replica::Pigeons.Replica{Vector{Float64}, NamedTuple{(:explorer_acceptance_pr, :explorer_n_steps, :swap_acceptance_pr, :log_sum_ratio, :traces, :index_process, :round_trip, :timing_extrema, :allocation_extrema, :_transformed_online), Tuple{OnlineStatsBase.GroupBy{Int64, Number, OnlineStatsBase.Mean{Float64, OnlineStatsBase.EqualWeight}}, OnlineStatsBase.GroupBy{Int64, Number, OnlineStatsBase.Sum{Float64}}, OnlineStatsBase.GroupBy{Tuple{Int64, Int64}, Number, OnlineStatsBase.Mean{Float64, OnlineStatsBase.EqualWeight}}, OnlineStatsBase.GroupBy{Tuple{Int64, Int64}, Number, Pigeons.LogSum{Float64}}, Dict{Pair{Int64, Int64}, Any}, Dict{Int64, Vector{Int64}}, Pigeons.RoundTripRecorder, NonReproducible{...}, NonReproducible{...}, Pigeons.OnlineStateRecorder}}}, pointer::Base.RefArray{Float64, Vector{Float64}, Nothing}, log_potential::InterpolatedLogPotential{...}, cached_lp::Float64)
    @ Pigeons ~/.julia/packages/Pigeons/rBo9q/src/explorers/SliceSampler.jl:81
  [3] slice_sample!(h::Pigeons.SliceSampler, state::Vector{Float64}, log_potential::InterpolatedLogPotential{...}, cached_lp::Float64, replica::Pigeons.Replica{Vector{Float64}, NamedTuple{(:explorer_acceptance_pr, :explorer_n_steps, :swap_acceptance_pr, :log_sum_ratio, :traces, :index_process, :round_trip, :timing_extrema, :allocation_extrema, :_transformed_online), Tuple{OnlineStatsBase.GroupBy{Int64, Number, OnlineStatsBase.Mean{Float64, OnlineStatsBase.EqualWeight}}, OnlineStatsBase.GroupBy{Int64, Number, OnlineStatsBase.Sum{Float64}}, OnlineStatsBase.GroupBy{Tuple{Int64, Int64}, Number, OnlineStatsBase.Mean{Float64, OnlineStatsBase.EqualWeight}}, OnlineStatsBase.GroupBy{Tuple{Int64, Int64}, Number, Pigeons.LogSum{Float64}}, Dict{Pair{Int64, Int64}, Any}, Dict{Int64, Vector{Int64}}, Pigeons.RoundTripRecorder, NonReproducible{...}, NonReproducible{...}, Pigeons.OnlineStateRecorder}}})
    @ Pigeons ~/.julia/packages/Pigeons/rBo9q/src/explorers/SliceSampler.jl:49
  [4] step!(explorer::Pigeons.SliceSampler, replica::Pigeons.Replica{Vector{Float64}, NamedTuple{(:explorer_acceptance_pr, :explorer_n_steps, :swap_acceptance_pr, :log_sum_ratio, :traces, :index_process, :round_trip, :timing_extrema, :allocation_extrema, :_transformed_online), Tuple{OnlineStatsBase.GroupBy{Int64, Number, OnlineStatsBase.Mean{Float64, OnlineStatsBase.EqualWeight}}, OnlineStatsBase.GroupBy{Int64, Number, OnlineStatsBase.Sum{Float64}}, OnlineStatsBase.GroupBy{Tuple{Int64, Int64}, Number, OnlineStatsBase.Mean{Float64, OnlineStatsBase.EqualWeight}}, OnlineStatsBase.GroupBy{Tuple{Int64, Int64}, Number, Pigeons.LogSum{Float64}}, Dict{Pair{Int64, Int64}, Any}, Dict{Int64, Vector{Int64}}, Pigeons.RoundTripRecorder, NonReproducible{...}, NonReproducible{...}, Pigeons.OnlineStateRecorder}}}, shared::Shared{...})
    @ Pigeons ~/.julia/packages/Pigeons/rBo9q/src/explorers/SliceSampler.jl:28
  [5] explore!(pt::PT{...}, replica::Pigeons.Replica{Vector{Float64}, NamedTuple{(:explorer_acceptance_pr, :explorer_n_steps, :swap_acceptance_pr, :log_sum_ratio, :traces, :index_process, :round_trip, :timing_extrema, :allocation_extrema, :_transformed_online), Tuple{OnlineStatsBase.GroupBy{Int64, Number, OnlineStatsBase.Mean{Float64, OnlineStatsBase.EqualWeight}}, OnlineStatsBase.GroupBy{Int64, Number, OnlineStatsBase.Sum{Float64}}, OnlineStatsBase.GroupBy{Tuple{Int64, Int64}, Number, OnlineStatsBase.Mean{Float64, OnlineStatsBase.EqualWeight}}, OnlineStatsBase.GroupBy{Tuple{Int64, Int64}, Number, Pigeons.LogSum{Float64}}, Dict{Pair{Int64, Int64}, Any}, Dict{Int64, Vector{Int64}}, Pigeons.RoundTripRecorder, NonReproducible{...}, NonReproducible{...}, Pigeons.OnlineStateRecorder}}}, explorer::Pigeons.SliceSampler)
    @ Pigeons ~/.julia/packages/Pigeons/rBo9q/src/pt/pigeons.jl:107
  [6] explore!
    @ ~/.julia/packages/Pigeons/rBo9q/src/pt/pigeons.jl:96 [inlined]
  [7] macro expansion
    @ ~/.julia/packages/Pigeons/rBo9q/src/pt/pigeons.jl:50 [inlined]
  [8] macro expansion
    @ ./timing.jl:501 [inlined]
  [9] run_one_round!(pt::PT{...})
    @ Pigeons ~/.julia/packages/Pigeons/rBo9q/src/pt/pigeons.jl:49
 [10] pigeons(pt::PT{...})
    @ Pigeons ~/.julia/packages/Pigeons/rBo9q/src/pt/pigeons.jl:18
 [11] pigeons(pt_arguments::Pigeons.Inputs{LogPosteriorDensity{Gravity.LikelihoodOptions{Float64, SourcePlane{Float64, Tuple{Vector{ImageFamily{Float64, PointImageMeasurement{Float64}}}}, Float64}, ImageFamily, ImageFamily{Float64, PointImageMeasurement{Float64}}, NamedTuple{(:bayesianfactor, :matches, :duplicates, :missedpredictionspenalty), Tuple{Bool, Symbol, Bool, Bool}}}}, Pigeons.GaussianReference, Pigeons.SliceSampler, LogPriorDensity{Gravity.LikelihoodOptions{Float64, SourcePlane{Float64, Tuple{Vector{ImageFamily{Float64, PointImageMeasurement{Float64}}}}, Float64}, ImageFamily, ImageFamily{Float64, PointImageMeasurement{Float64}}, NamedTuple{(:bayesianfactor, :matches, :duplicates, :missedpredictionspenalty), Tuple{Bool, Symbol, Bool, Bool}}}}}, #unused#::Pigeons.ThisProcess)
    @ Pigeons ~/.julia/packages/Pigeons/rBo9q/src/api.jl:19
 [12] #pigeons#157
    @ ~/.julia/packages/Pigeons/rBo9q/src/api.jl:16 [inlined]

Is there something I should check to avoid this error? From my tests, my logposterior seems to be finite and I do not see where the problem could be. Thank you for any help you can provide!

[alexandrebouchard]

This is often a symptom that the posterior is improper (due to an intentionally improper prior, or even I had it when forgetting to specify a prior on a variable). In such case the slice sampler might venture off very far where numerical instability could end up giving a NaN even when the target is well behaved.

If not, feel free to send instructions to replicate the issue.

[astrozot]

Thank you so much Alexandre for your suggestions! I will double and check if I can solve this...

[miguelbiron]

Hi @astrozot -- I have a PR (#183) with some improvements on the SliceSampler code, and in particular with a more detailed error whenever invalid log densities are encountered. If you want to test it, you can do this by adding this branch of Pigeons via

using Pkg
Pkg.add(name="Pigeons", rev="ss-handle-inf-z")

That should hopefully raise the exact point in the state space that is found to have invalid density.

[miguelbiron]

Just merged the PR so you can simply install the latest dev version, e.g. via

using Pkg
Pkg.add(name="Pigeons", rev="main")

[astrozot]

Hi @miguelbiron, thank you so much for your message! I will install the latest dev version and check if I can get a clearer clue of my issue.