Strange error from Turing.jl?

Specific Domains Probabilistic programming
1

I’ve got the following model, and some data, I’m trying to get a posterior over the infection probability for vaxxed and unvaxxed in my simulated dataset:


@model function vaxxdecision(data)

    pv ~ Beta(1,100)
    punv ~ Beta(1,100)

    vaxxed = sum(data.vaxxed==true)
    unvaxxed = sum(data.vaxxed==false)

    vaxxedinf = sum(data[data.vaxxed .==true,:infected])
    unvaxxedinf = sum(data[data.vaxxed .==false,:infected])

    vaxxedinf ~ Binomial(vaxxed,pv)
    unvaxxedinf ~ Binomial(unvaxxed,punv)
end


simdata = DataFrame(vaxxed = [s.vaxdose > 0 for s in sims],
            infected=[x.severity > .01 for x in sims])


chain = sample(vaxxdecision(simdata),NUTS(500,.85),1000)

The last line throws an error
ERROR: DomainError with (-0.23237667153547514, 2.232376671535475):
beta(a, b) must be non-negative

Long Stack Trace Elided 
ERROR: DomainError with (-0.23237667153547514, 2.232376671535475):
`beta(a, b)` must be non-negative
Stacktrace:
  [1] logbeta
    @ ~/.julia/packages/SpecialFunctions/tqwrL/src/gamma.jl:790 [inlined]
  [2] logbeta
    @ ~/.julia/packages/ForwardDiff/5gUap/src/dual.jl:233 [inlined]
  [3] binomlogpdf(n::Int64, p::ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#7"{DynamicPPL.TypedVarInfo{NamedTuple{(:pv, :punv, :vaxxedinf, :unvaxxedinf), Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:pv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:pv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:punv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:punv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:vaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:vaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{var"#87#88", (:data,), (), (), Tuple{DataFrame}, Tuple{}}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 4}, k::ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#7"{DynamicPPL.TypedVarInfo{NamedTuple{(:pv, :punv, :vaxxedinf, :unvaxxedinf), Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:pv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:pv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:punv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:punv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:vaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:vaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{var"#87#88", (:data,), (), (), Tuple{DataFrame}, Tuple{}}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 4})
    @ StatsFuns ~/.julia/packages/StatsFuns/vyFzQ/src/distrs/binom.jl:19
  [4] logpdf
    @ ~/.julia/packages/Distributions/XCZG8/src/univariates.jl:653 [inlined]
  [5] #5
    @ ~/.julia/packages/Bijectors/OJrCc/src/Bijectors.jl:151 [inlined]
  [6] map
    @ ./number.jl:272 [inlined]
  [7] logpdf_with_trans
    @ ~/.julia/packages/Bijectors/OJrCc/src/Bijectors.jl:151 [inlined]
  [8] assume(rng::Random._GLOBAL_RNG, 

... CUT HERE because the total stack trace is too long to be allowed for posting

Anyone have any idea what this error means?

2

Weird, I noticed an error where I wasn’t broadcasting when doing the counts of vaxxed and unvaxxed, and fixed it, now I have a different error:



@model function vaxxdecision(data)

    pv ~ Beta(1.0,100.0)
    punv ~ Beta(1.0,100.0)

    vaxxed = sum(data.vaxxed .==true)
    unvaxxed = sum(data.vaxxed .==false)

    vaxxedinf = sum(data[data.vaxxed .==true,:infected])
    unvaxxedinf = sum(data[data.vaxxed .==false,:infected])

    vaxxedinf ~ Binomial(vaxxed,pv)
    unvaxxedinf ~ Binomial(unvaxxed,punv)
end


simdata = DataFrame(vaxxed = [s.vaxdose > 0 for s in sims],
            infected=[x.severity > .01 for x in sims])




chain = sample(vaxxdecision(simdata),NUTS(500,.01),1000)

julia> chain = sample(vaxxdecision(simdata),NUTS(500,.85),1000)

┌ Info: Found initial step size
└   ϵ = 0.2
Sampling 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:02
ERROR: InexactError: Int64(359.01732013963164)
Stacktrace:
  [1] Int64
    @ ./float.jl:812 [inlined]
  [2] convert
    @ ./number.jl:7 [inlined]
  [3] setindex!
    @ ./array.jl:877 [inlined]
  [4] setindex!
    @ ./subarray.jl:317 [inlined]
  [5] copyto_unaliased!(deststyle::IndexCartesian, dest::SubArray{Int64, 1, Vector{Int64}, Tuple{Vector{Int64}}, false}, srcstyle::IndexLinear, src::SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true})
    @ Base ./abstractarray.jl:1032
  [6] copyto!
    @ ./abstractarray.jl:998 [inlined]
  [7] copyto!
    @ ./broadcast.jl:998 [inlined]
  [8] copyto!
    @ ./broadcast.jl:957 [inlined]
  [9] materialize!
    @ ./broadcast.jl:915 [inlined]
 [10] materialize!
    @ ./broadcast.jl:912 [inlined]
 [11] macro expansion
    @ ~/.julia/packages/DynamicPPL/wf0dU/src/varinfo.jl:0 [inlined]
 [12] _setindex!(metadata::NamedTuple{(:pv, :punv, :vaxxedinf, :unvaxxedinf), Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:pv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:pv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:punv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:punv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:vaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:vaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}}}, val::Vector{Float64}, ranges::NamedTuple{(:pv, :punv, :vaxxedinf, :unvaxxedinf), NTuple{4, Vector{Int64}}})
    @ DynamicPPL ~/.julia/packages/DynamicPPL/wf0dU/src/varinfo.jl:892
 [13] setindex!(vi::DynamicPPL.TypedVarInfo{NamedTuple{(:pv, :punv, :vaxxedinf, :unvaxxedinf), Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:pv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:pv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:punv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:punv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:vaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:vaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, val::Vector{Float64}, spl::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}})
    @ DynamicPPL ~/.julia/packages/DynamicPPL/wf0dU/src/varinfo.jl:888
 [14] step(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{var"#115#116", (:data,), (), (), Tuple{DataFrame}, Tuple{}}, spl::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, state::Turing.Inference.HMCState{DynamicPPL.TypedVarInfo{NamedTuple{(:pv, :punv, :vaxxedinf, :unvaxxedinf), Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:pv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:pv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:punv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:punv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:vaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:vaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, AdvancedHMC.NUTS{AdvancedHMC.MultinomialTS, AdvancedHMC.GeneralisedNoUTurn, AdvancedHMC.Leapfrog{Float64}, Float64}, AdvancedHMC.Hamiltonian{AdvancedHMC.DiagEuclideanMetric{Float64, Vector{Float64}}, Turing.Inference.var"#logπ#54"{DynamicPPL.TypedVarInfo{NamedTuple{(:pv, :punv, :vaxxedinf, :unvaxxedinf), Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:pv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:pv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:punv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:punv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:vaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:vaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.Model{var"#115#116", (:data,), (), (), Tuple{DataFrame}, Tuple{}}}, Turing.Inference.var"#∂logπ∂θ#53"{DynamicPPL.TypedVarInfo{NamedTuple{(:pv, :punv, :vaxxedinf, :unvaxxedinf), Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:pv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:pv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:punv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:punv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:vaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:vaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.Model{var"#115#116", (:data,), (), (), Tuple{DataFrame}, Tuple{}}}}, AdvancedHMC.PhasePoint{Vector{Float64}, AdvancedHMC.DualValue{Float64, Vector{Float64}}}, AdvancedHMC.Adaptation.StanHMCAdaptor{AdvancedHMC.Adaptation.WelfordVar{Float64, Vector{Float64}}, AdvancedHMC.Adaptation.NesterovDualAveraging{Float64}}}; nadapts::Int64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Turing.Inference ~/.julia/packages/Turing/uAz5c/src/inference/hmc.jl:292
 [15] macro expansion
    @ ~/.julia/packages/AbstractMCMC/oou1a/src/sample.jl:108 [inlined]
 [16] macro expansion
    @ ~/.julia/packages/ProgressLogging/6KXlp/src/ProgressLogging.jl:328 [inlined]
 [17] (::AbstractMCMC.var"#20#21"{Bool, String, Nothing, Int64, Int64, Base.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:nadapts,), Tuple{Int64}}}, Random._GLOBAL_RNG, DynamicPPL.Model{var"#115#116", (:data,), (), (), Tuple{DataFrame}, Tuple{}}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, Int64, Int64})()
    @ AbstractMCMC ~/.julia/packages/AbstractMCMC/oou1a/src/logging.jl:11
 [18] with_logstate(f::Function, logstate::Any)
...
3

Tried rewriting it like this, but still get similar error:


@model function vaxxdecision(data)

    pv ~ Beta(1.0,100.0)
    punv ~ Beta(1.0,100.0)

    vaxxed = 0
    unvaxxed=0
    vaxxedinf = 0
    unvaxxedinf = 0

    for i in 1:nrow(data)
        if data[i,:vaxxed] == true
            vaxxed +=1
            if data[i,:infected] == true
                vaxxedinf +=1
            end 
        else
            unvaxxed += 1
            if data[i,:infected] == true
                unvaxxedinf +=1
            end 

        end            
    end

    vaxxedinf ~ Binomial(vaxxed,pv)
    unvaxxedinf ~ Binomial(unvaxxed,punv)
end


simdata = DataFrame(vaxxed = [s.vaxdose > 0 for s in sims],
            infected=[x.severity > .01 for x in sims])




chain = sample(vaxxdecision(simdata),NUTS(500,.85;init_ϵ=.001),1000)

error:

ERROR: InexactError: Int64(72.89246517446986)
Stacktrace:
  [1] Int64
    @ ./float.jl:812 [inlined]
  [2] convert
    @ ./number.jl:7 [inlined]
  [3] setindex!
    @ ./array.jl:877 [inlined]
  [4] setindex!
    @ ./subarray.jl:317 [inlined]
  [5] copyto_unaliased!(deststyle::IndexCartesian, dest::SubArray{Int64, 1, Vector{Int64}, Tuple{Vector{Int64}}, false}, srcstyle::IndexLinear, src::SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true})
    @ Base ./abstractarray.jl:1032
  [6] copyto!
    @ ./abstractarray.jl:998 [inlined]
  [7] copyto!
    @ ./broadcast.jl:998 [inlined]
  [8] copyto!
    @ ./broadcast.jl:957 [inlined]
  [9] materialize!
    @ ./broadcast.jl:915 [inlined]
 [10] materialize!
    @ ./broadcast.jl:912 [inlined]
 [11] macro expansion
    @ ~/.julia/packages/DynamicPPL/wf0dU/src/varinfo.jl:0 [inlined]
 [12] _setindex!(metadata::NamedTuple{(:pv, :punv, :vaxxedinf, :unvaxxedinf), Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:pv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:pv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:punv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:punv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:vaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:vaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}}}, val::Vector{Float64}, ranges::NamedTuple{(:pv, :punv, :vaxxedinf, :unvaxxedinf), NTuple{4, Vector{Int64}}})
    @ DynamicPPL ~/.julia/packages/DynamicPPL/wf0dU/src/varinfo.jl:892
 [13] setindex!(vi::DynamicPPL.TypedVarInfo{NamedTuple{(:pv, :punv, :vaxxedinf, :unvaxxedinf), Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:pv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:pv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:punv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:punv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:vaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:vaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, val::Vector{Float64}, spl::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}})
    @ DynamicPPL ~/.julia/packages/DynamicPPL/wf0dU/src/varinfo.jl:888
 [14] initialstep(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{var"#133#134", (:data,), (), (), Tuple{DataFrame}, Tuple{}}, spl::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, vi::DynamicPPL.TypedVarInfo{NamedTuple{(:pv, :punv, :vaxxedinf, :unvaxxedinf), Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:pv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:pv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:punv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:punv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:vaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:vaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}; init_params::Nothing, nadapts::Int64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Turing.Inference ~/.julia/packages/Turing/uAz5c/src/inference/hmc.jl:218
 [15] step(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{var"#133#134", (:data,), (), (), Tuple{DataFrame}, Tuple{}}, spl::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}; resume_from::Nothing, kwargs::Base.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:nadapts,), Tuple{Int64}}})
4

I have no idea why it is doing that, but I get the same error when I try use simulated data. You could try pulling out the data and summarizing it before putting it into the model, like in the code example below:

Random.seed!(98245)
sims = 1000
p_vax = 0.75
p_inf_vax = 0.2
p_inf_unvax = 0.8
vaxxed = rand(sims) .< p_vax
infected = zeros(Int64, sims)

for i in 1:sims
    p = rand() < p_inf_vax 
    t = vaxxed[i] == 1 
    
    if t == true & p == true
        infected[i] = 1
    end
end


for i in 1:sims
    p = rand() < p_inf_unvax 
    t = vaxxed[i] == 0

    if t == true & p == true
        infected[i] = 1
    end

end
        

simdata = DataFrame(:vaxxed => [vaxxed[i] .== 1 for i in 1:sims],
                    :infected => [infected[i] .== 1 for i in 1:sims])


vaxxed = sum(simdata.vaxxed .== true)
unvaxxed = sum(simdata.vaxxed .== false)
vaxxedinf = sum(simdata.infected[simdata.vaxxed .== true])
unvaxxedinf = sum(simdata.infected[simdata.vaxxed .== false])


@model function vaxxdecision4(vaxxed, unvaxxed, vaxxedinf, unvaxxedinf)

    pv ~ Beta(1,100)
    punv ~ Beta(1,100)

    vaxxedinf ~ Binomial(vaxxed, pv)
    unvaxxedinf ~ Binomial(unvaxxed, punv)

end

chn = sample(vaxxdecision4(vaxxed, unvaxxed, vaxxedinf, unvaxxedinf), NUTS(), 1000)
5

yeah, the issue is that in the “real” model, the plan is to have more complicated stuff. If i can’t do this very simple thing, what the heck can I do?

6

In general, you want to do all of your data processing outside of the Turing model, since any done inside the model will need to be repeated every time the log density is evaluated. So try something like this:

@model function vaxxdecision(vaxxed, unvaxxed, vaxxedinf, unvaxxedinf)
    pv ~ Beta(1,100)
    punv ~ Beta(1,100)

    vaxxedinf ~ Binomial(vaxxed,pv)
    unvaxxedinf ~ Binomial(unvaxxed,punv)
end

vaxxed = sum(data.vaxxed==true)
unvaxxed = sum(data.vaxxed==false)

vaxxedinf = sum(data[data.vaxxed .==true,:infected])
unvaxxedinf = sum(data[data.vaxxed .==false,:infected])

chain = sample(
    vaxxdecision(vaxxed, unvaxxed, vaxxedinf, unvaxxedinf),
    NUTS(500, .85),
    1000,
)

As the provided example is not complete, I cannot verify that this works. But in this particular case, the posterior is known to be a Beta distribution.

Specifically, the posterior should be something like

pv ~ Beta(1 + vaxxedinf, 100 + vaxxed - vaxxedinf)
punv ~ Beta(1 + unvaxxedinf, 100 + unvaxxed - unvaxxedinf)
1 Like
7

Ok, makes sense. Yes, this is just a toy model, the real model will be significantly more complicated. It’s intended to be a decision theoretic model using frequency of various levels of side effects, and frequency of various main effects of the disease, and various other things. It’ll condition on various things for each patient, so I don’t think I can pre-compute everything for the real model as they aren’t all homogeneous.

For example, side effects could be a function of age, weight, height, and dose administered. all of which are unique to each patient.

8

What’s interesting to me is that when I tried to duplicate your original model, I got the same error:

julia> chn = sample(vaxxdecision(simdata), NUTS(), 1000)
┌ Info: Found initial step size
└   ϵ = 0.2
ERROR: InexactError: Int64(608.9873768061884)
Stacktrace:
  [1] Int64

But in the simulated data I created, none of those numbers should add up to ~609 (in this case, anyway. I get a different number each time I try rerun it. So far the range had been 314 to 609. But my point remains the same).

julia> vaxxed = sum(simdata.vaxxed .== true)
722

julia> unvaxxed = sum(simdata.vaxxed .== false)
278

julia> vaxxedinf = sum(simdata.infected[simdata.vaxxed .== true])
142

julia> unvaxxedinf = sum(simdata.infected[simdata.vaxxed .== false])
214

I have no idea where the inexact sums are coming from.

On a separate note, I also tried rewriting the model to be type-stable, but still got similar errors.

@model function vaxxdecision3(data, ::Type{T} = Int64) where {T} 

    pv ~ Beta(1, 100)
    punv ~ Beta(1, 100)

    numinf = Vector{T}(undef, 2)
    vaxxed = Vector{T}(undef, 2)

    numinf[1] = sum(data[data.vaxxed .== true, :infected])
    numinf[2] = sum(data[data.vaxxed .== false, :infected])

    vaxxed[1] = sum(data.vaxxed .== true)
    vaxxed[2] = sum(data.vaxxed .== false)

    numinf[1] ~ Binomial(vaxxed[1], pv)
    numinf[2] ~ Binomial(vaxxed[2], punv)
end

chn = sample(vaxxdecision3(simdata), NUTS(), 1000)

ERROR: InexactError: Int(Int64, Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:pv, :punv, :numinf), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:pv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{AbstractPPL.VarName{:pv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:punv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{AbstractPPL.VarName{:punv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:numinf, Tuple{Tuple{Int64}}}, Int64}, Vector{Binomial{Float64}}, Vector{AbstractPPL.VarName{:numinf, Tuple{Tuple{Int64}}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(vaxxdecision3), (:data, :T), (:T,), (), Tuple{DataFrame, DataType}, Tuple{DataType}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}}(437.0,0.0,0.0,1.0,0.0))
Stacktrace:
  [1] Int64
...

I leave this to someone with a better understanding of Turing, DynamicPPL, etc., to try figure out.

9

I haven’t looked at the code indicated by my stacktrace yet but it was within a setindex! Call I wonder if this is a bug where it’s trying to index something with a float?

10

Have you tried the suggestions in this comment yet? Did they work? If not, please share the stacktrace.

11

Unfortunately that approach, even if it solves the problem doesn’t really solve the real problem, since the real model will not use a simple summary but will have to loop through all the observations.

@sethaxen
I tried the following and it works, so my Turing install isn’t totally broken. Unfortunately as I said above, the real model isn’t like this, it’s more a regression where we predict outcomes for each patient based on their individual characteristics.


@model function simplevax(nv,nvinf,nnv,nnvinf)
    pv ~ Beta(1.0,100.0)
    punv ~ Beta(1.0,100.0)
    nvinf ~ Binomial(nv,pv)
    nnvinf ~ Binomial(nnv,punv)
end

chaintest = sample(simplevax(500,3,500,30),NUTS(500,.85;init_ϵ=.001),1000)
12

I see. Can you share a complete minimal failing example? In your original post sims is not defined, so we don’t have all of the information to debug with.

13

Sure, sims is created from the output of a simulation:
But for purposes of reproducing the issue minimally you could just do this for simdata:

simdata = DataFrame(vaxxed = [ i < 500 for i from 1:1000], infected = rand(Bernoulli(.02),1000))

When I look at the backtrace, eventually I get here…

 [13] setindex!(vi::DynamicPPL.TypedVarInfo{NamedTuple{(:pv, :punv, :vaxxedinf, :unvaxxedinf), Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:pv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:pv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:punv, Tuple{}}, Int64}, Vector{Beta{Float64}}, Vector{DynamicPPL.VarName{:punv, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:vaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:vaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}, Int64}, Vector{Binomial{Float64}}, Vector{DynamicPPL.VarName{:unvaxxedinf, Tuple{}}}, Vector{Int64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, val::Vector{Float64}, spl::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}})

I have no idea what this really means, but I wonder if:

DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:pv, Tuple{}}, Int64},

suggests that it thinks pv is a variable with type Int64, even though I’m clearly giving it a prior with a Beta distribution?

If it’s helpful to have it all in one place, here’s a MWE:


@model function vaxxdecision(data)

    pv ~ Beta(1.0,100.0)
    punv ~ Beta(1.0,100.0)

    vaxxed = 0
    unvaxxed=0
    vaxxedinf = 0
    unvaxxedinf = 0

    for i in 1:nrow(data)
        if data[i,:vaxxed] == true
            vaxxed +=1
            if data[i,:infected] == true
                vaxxedinf +=1
            end 
        else
            unvaxxed += 1
            if data[i,:infected] == true
                unvaxxedinf +=1
            end 

        end            
    end

    vaxxedinf ~ Binomial(vaxxed,pv)
    unvaxxedinf ~ Binomial(unvaxxed,punv)
end


simdata = DataFrame(vaxxed = [i < 500 for i in 1:1000],infected = rand(Bernoulli(.02),1000))


chain = sample(vaxxdecision(simdata),NUTS(500,.85;init_ϵ=.001),1000)
14

Ok, I thought maybe it was something specific to the NUTS sampler, so I tried with SMC() just for giggles, and I got a Segfault:

julia> chain = sample(vaxxdecision(simdata),SMC(),1000)

signal (11): Segmentation fault
in expression starting at REPL[3]:1
setproperty! at ./Base.jl:43 [inlined]
consume at /home/dlakelan/.julia/packages/Libtask/OtgFl/src/ctask.jl:180
consume at /home/dlakelan/.julia/packages/Turing/uAz5c/src/core/container.jl:57 [inlined]
reweight! at /home/dlakelan/.julia/packages/Turing/uAz5c/src/core/container.jl:252
sweep! at /home/dlakelan/.julia/packages/Turing/uAz5c/src/core/container.jl:307
#initialstep#73 at /home/dlakelan/.julia/packages/Turing/uAz5c/src/inference/AdvancedSMC.jl:121
initialstep##kw at /home/dlakelan/.julia/packages/Turing/uAz5c/src/inference/AdvancedSMC.jl:111
unknown function (ip: 0x7fa4a449a08f)
_jl_invoke at /buildworker/worker/package_linux64/build/src/gf.c:2245 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:2427
#step#18 at /home/dlakelan/.julia/packages/DynamicPPL/wf0dU/src/sampler.jl:83
step##kw at /home/dlakelan/.julia/packages/DynamicPPL/wf0dU/src/sampler.jl:66 [inlined]
macro expansion at /home/dlakelan/.julia/packages/AbstractMCMC/oou1a/src/sample.jl:97 [inlined]
macro expansion at /home/dlakelan/.julia/packages/ProgressLogging/6KXlp/src/ProgressLogging.jl:328 [inlined]
#20 at /home/dlakelan/.julia/packages/AbstractMCMC/oou1a/src/logging.jl:11
unknown function (ip: 0x7fa4a448eeef)
_jl_invoke at /buildworker/worker/package_linux64/build/src/gf.c:2245 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:2427
with_logstate at ./logging.jl:511
with_logger at ./logging.jl:623
unknown function (ip: 0x7fa4a496e681)
_jl_invoke at /buildworker/worker/package_linux64/build/src/gf.c:2245 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:2427
with_progresslogger at /home/dlakelan/.julia/packages/AbstractMCMC/oou1a/src/logging.jl:34
unknown function (ip: 0x7fa4a496dc98)
_jl_invoke at /buildworker/worker/package_linux64/build/src/gf.c:2245 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:2427
macro expansion at /home/dlakelan/.julia/packages/AbstractMCMC/oou1a/src/logging.jl:10 [inlined]
#mcmcsample#19 at /home/dlakelan/.julia/packages/AbstractMCMC/oou1a/src/sample.jl:88
unknown function (ip: 0x7fa4a4482fee)
mcmcsample##kw at /home/dlakelan/.julia/packages/AbstractMCMC/oou1a/src/sample.jl:85

However if I try the metropolis hastings algorithm MH(), it works:

chain = sample(vaxxdecision(simdata),MH(),1000)

Strangely though, it treats vaxxedinf and unvaxxedinf as if they were parameters, rather than data:

image
image915×951 149 KB

15

AHA! I have figured it out by doing some @macroexpand calls on the @model macro.

Turing looks at the arguments in the model function to determine whether something on the left hand side of ~ is an observation or something it should generate a sample for…

So this works:


@model function vaxxdecision(data,vaxxedinf, unvaxxedinf)

    pv ~ Beta(1.0,100.0)
    punv ~ Beta(1.0,100.0)

    vaxxed = 0
    unvaxxed=0
    vaxxedinf = 0
    unvaxxedinf = 0

    for i in 1:nrow(data)
        if data[i,:vaxxed] == true
            vaxxed +=1
            if data[i,:infected] == true
                vaxxedinf +=1
            end 
        else
            unvaxxed += 1
            if data[i,:infected] == true
                unvaxxedinf +=1
            end 

        end            
    end

    vaxxedinf ~ Binomial(vaxxed,pv)
    unvaxxedinf ~ Binomial(unvaxxed,punv)
end

simdata = DataFrame(vaxxed = [s.vaxdose > 0 for s in sims],
            infected=[x.severity > .01 for x in sims])

chain = sample(vaxxdecision(simdata,0,0),NUTS(500,.85;init_ϵ=.001),1000)

And it doesn’t sample from the vaxxedinf and unvaxxedinf, it just uses pv and punv as the two parameters:

image
image939×522 85.7 KB

1 Like
16

Ah, nice work!

2 Likes
17

It would be great to have this very explicitly somewhere early-on in the Turing documentation.

Anyone have suggestions as to how we can make this super clear in the docs?

18

this seems to be the same error I get when using SMC. This can be reproduced on my machine even with the gdemo example from the docs. I created an issue: Turing segfaults on gdemo example with SMC on julia · Issue #1737 · TuringLang/Turing.jl · GitHub

19

I guess based on above discussion that the segfault is expected on 1.7 due to libtask issues (used in all of the “particle” methods)