# Turing model fail due to extreme outliers sampled in Exponential distribution?

I am modelling some textbook models (minimal example given below)

using Turing

# Data ===================
contact = [1, 1, 1, 2, 2, 2, 2, 1, 2, 1]
total_tools = [ 13 ,22 ,24 ,43 ,33 ,19 ,40 ,28 ,55 ,71]
population = [1100, 1500, 3600, 4791, 7400, 8000, 9200, 13000, 17500, 275000]
population = (population .- mean(population))./std(population) # normalize
population = population .+ abs(minimum(population)) .+ 0.1    # shift to positive values

# model =============================
@model function m12_2(tools, population, cid)
ϕ ~ Exponential(1.0)
g ~ Exponential(1.0)
a = Vector{Real}(undef,2)
b = Vector{Real}(undef,2)
a .~ Normal(1.0, 1.0)
b .~ Exponential(1.0)
for i = 1 : length(cid)
λ = (exp(a[cid[i]]) * population[i]^ b[cid[i]])/g
tools[i] ~ NegativeBinomial(λ/(ϕ + 0.000001), 1.0 / (1.0 + ϕ))
end
end
describe(sample(m12_2(total_tools, population, contact), NUTS(),1000))


Problem is that above model fails randomly 2 out 3 times with error message:

ArgumentError: NegativeBinomial: the condition r > zero(r) is not satisfied.

macro expansion@utils.jl:6[inlined]
#NegativeBinomial#46@negativebinomial.jl:42[inlined]
NegativeBinomial@negativebinomial.jl:41[inlined]
#1@Other: 10[inlined]
(::Main.workspace#160.var"#1#2")(::DynamicPPL.Model{Main.workspace#160.var"#1#2", (:tools, :population, :cid), (), (), Tuple{Vector{Int64}, Vector{Float64}, Vector{Int64}}, Tuple{}}, ::DynamicPPL.ThreadSafeVarInfo{DynamicPPL.TypedVarInfo{NamedTuple{(:ϕ, :g, :a, :b), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ϕ, Tuple{}}, Int64}, Vector{Distributions.Exponential{Float64}}, Vector{AbstractPPL.VarName{:ϕ, Tuple{}}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#3"{DynamicPPL.TypedVarInfo{NamedTuple{(:ϕ, :g, :a, :b), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ϕ, Tuple{}}, Int64}, Vector{Distributions.Exponential{Float64}}, Vector{AbstractPPL.VarName{:ϕ, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:g, Tuple{}}, Int64}, Vector{Distributions.Exponential{Float64}}, Vector{AbstractPPL.VarName{:g, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, ...
...
Model@model.jl:91[inlined]
Model@model.jl:104[inlined]
vector_mode_dual_eval@apiutils.jl:37[inlined]
gradient_logp(::Turing.Core.ForwardDiffAD{40}, ::Vector{Float64}, ::DynamicPPL.TypedVarInfo{NamedTuple{(:ϕ, :g, :a, :b), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ϕ, Tuple{}}, Int64}, Vector{Distributions.Exponential{Float64}}, Vector{AbstractPPL.VarName{:ϕ, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:g, Tuple{}}, Int64}, Vector{Distributions.Exponential{Float64}}, Vector{AbstractPPL.VarName{:g, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:a, Tuple{Tuple{Int64}}}, Int64}, Vector{Distributions.Normal{Float64}}, Vector{AbstractPPL.VarName{:a, Tuple{Tuple{Int64}}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:b, Tuple{Tuple{Int64}}}, Int64}, Vector{Distributions.Exponential{Float64}}, Vector{AbstractPPL.VarName{:b, Tuple{Tuple{Int64}}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, ::DynamicPPL.Model{Main.workspace#160.var"#1#2", (:tools, :population, :cid), (), (), Tuple{Vector{Int64}, Vector{Float64}, Vector{Int64}}, Tuple{}}, ::DynamicPPL.Sampler{Turing.Inference.NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, ::DynamicPPL.DefaultContext)@ad.jl:121
∂logπ∂θ@hmc.jl:433[inlined]
∂H∂θ@hamiltonian.jl:31[inlined]
macro expansion@UnPack.jl:100[inlined]
step@integrator.jl:66[inlined]
transition@sampler.jl:57[inlined]
macro expansion@sample.jl:134[inlined]
macro expansion@ProgressLogging.jl:328[inlined]
macro expansion@logging.jl:8[inlined]
var"#mcmcsample#20"(::Bool, ::String, ::Nothing, ::Int64, ::Int64, ::Type, ::Base.Iterators.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:nadapts,), Tuple{Int64}}}, ::typeof(AbstractMCMC.mcmcsample), ::Random._GLOBAL_RNG, ::DynamicPPL.Model{Main.workspace#160.var"#1#2", (:tools, :population, :cid), (), (), Tuple{Vector{Int64}, Vector{Float64}, Vector{Int64}}, Tuple{}}, ::DynamicPPL.Sampler{Turing.Inference.NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, ::Int64)@sample.jl:114
var"#sample#40"(::Type, ::Nothing, ::Bool, ::Int64, ::Bool, ::Int64, ::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, ::typeof(StatsBase.sample), ::Random._GLOBAL_RNG, ::DynamicPPL.Model{Main.workspace#160.var"#1#2", (:tools, :population, :cid), (), (), Tuple{Vector{Int64}, Vector{Float64}, Vector{Int64}}, Tuple{}}, ::DynamicPPL.Sampler{Turing.Inference.NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, ::Int64)@hmc.jl:133
sample@hmc.jl:116[inlined]
#sample#2@Inference.jl:142[inlined]
sample@Inference.jl:142[inlined]
#sample#1@Inference.jl:132[inlined]
sample(::DynamicPPL.Model{Main.workspace#160.var"#1#2", (:tools, :population, :cid), (), (), Tuple{Vector{Int64}, Vector{Float64}, Vector{Int64}}, Tuple{}}, ::Turing.Inference.NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}, ::Int64)@Inference.jl:132
top-level scope@Local: 1[inlined]


To figure out the bug I used try ... catch  blocks inside my model as

@model function m12_2(tools, population, cid)
ϕ ~ Exponential(1.0)
g ~ Exponential(1.0)
a = Vector{Real}(undef,2)
b = Vector{Real}(undef,2)
a .~ Normal(1.0, 1.0)
b .~ Exponential(1.0)
for i = 1 : length(cid)
λ = (exp(a[cid[i]]) * population[i]^ b[cid[i]])/g
try
tools[i] ~ NegativeBinomial(λ/(ϕ + 0.000001), 1.0 / (1.0 + ϕ))
catch
println(i)
println(population[i], " ",b[cid[i]])
return -1
end
end
end


It looks like the random variable b, drawn from Exponential(1.0) sample extremely large values for some reason (like 2.0947637423773094e13 ?), resulting in population[i]^b[cid] becoming numerically zero.
Given below is the output of above model.



From worker 3:	1
From worker 3:	0.1 Dual{ForwardDiff.Tag{Turing.Core.var"#f#3"{DynamicPPL.TypedVarInfo{NamedTuple{(:ϕ, :g, :a, :b), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ϕ, Tuple{}}, Int64}, Vector{Distributions.Exponential{Float64}}, Vector{AbstractPPL.VarName{:ϕ, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:g, Tuple{}}, Int64}, Vector{Distributions.Exponential{Float64}}, Vector{AbstractPPL.VarName{:g, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:a, Tuple{Tuple{Int64}}}, Int64}, Vector{Distributions.Normal{Float64}}, Vector{AbstractPPL.VarName{:a, Tuple{Tuple{Int64}}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:b, Tuple{Tuple{Int64}}}, Int64}, Vector{Distributions.Exponential{Float64}}, Vector{AbstractPPL.VarName{:b, Tuple{Tuple{Int64}}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{Main.workspace#173.var"#1#2", (:tools, :population, :cid), (), (), Tuple{Vector{Int64}, Vector{Float64}, Vector{Int64}}, Tuple{}}, DynamicPPL.Sampler{Turing.Inference.NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}}(2.0947637423773094e13,0.0,0.0,0.0,0.0,2.0947637423773094e13,0.0)

1. Is it some bug or something wrong in my model?
2. Is there some sort of Turing debugger to assist in such matter?

Although there are plans to remove it eventually for now I think you should use filldist or arraydist to declare a and b. See the Turing docs

Is expected outcome of filldist different from .~ ? I tried it but the problem persist. Currently my only workaround is to use truncated(Exponential(1.0), 1, 10) .

Which sampling algorithm are you using?

NUTS. Though I tried SMC, and HMC as well. But the same error persisted.

So this is sampling as a parameter, if it’s heading off to infinity it might be because some other parameter collapses and you get divergent transitions. try NUTS(1000,0.99) to make it do a very careful solution to the HMC differential equations, this might help?

you might also try making this fudge factor larger, like phi + 0.001

I made fudge factor = 0.001 and sampled using NUTS(1000,0.99) but the problem still persists.

weird. does this occur during the NUTS warmup or after?

I would assume its after warmup, as it only gives the error if i draw samples > 500 or so. At lower samples it rarely crashes if at all.

I suspect if phi goes to zero then lambda needs to go to zero to maintain a reasonable ratio. Try making the fudge factor something like 0.2 just to see if it fixes the issue.

No. Still the same with fudge factor of 0.2.
Only truncated(Exponential()) works as of now.

Does it do the same thing with RWMH algorithm?

Ok, so just to make sure there isn’t a basic problem with sampling I ran this:

using Turing, StatsPlots, MCMCChains

@model function f()
x ~ Exponential(1.0)
end

ch = sample(f(),NUTS(1000,.8),2000)


And it sampled normally with 97.5%tile of 3.8621 under latest version of Turing

Then I copy pasta your original code and get the error

ERROR: LoadError: ArgumentError: NegativeBinomial: the condition r > zero(r) is not satisfied.
Stacktrace:


The fact that this is an intermittent bug suggests to me that the model itself is somehow unstable.

If I modify the code as follows:

		λ = (exp(a[cid[i]]) * (.00001 + population[i]^ b[cid[i]]))/g


The model samples just fine and doesn’t give strange values for b. I think this is a modeling issue rather than a sampling issue

When you would modify it as (.00001 + population[i]^ b[cid[i]])) issue will be mitigated as \lambda will no longer go to zero. Problem is that population[i]^ b[cid[i]] goes to zero if b goes higher. With addition of 0.0001, values of b will no longer reduce \lambda to zero.

To replicate it in simpler model i tried the following:

@model function f()
x ~ Exponential(1.0)
y = 0.1^x
z ~ NegativeBinomial(1/y,0.1)
end


But doing so does not reproduce the original error. I do get warnings though.

┌ Warning: The current proposal will be rejected due to numerical error(s).
│   isfinite.((θ, r, ℓπ, ℓκ)) = (true, true, false, true)