Making Turing Fast with large numbers of parameters?

General Usage
81

I tried this:


function truncatenormal(a,b)::UnivariateDistribution
    truncated(Normal(a,b),0.0,Inf)
end


@model function estweights3lazy(nstaff,staffid,npatients,patientid,usewch,totweight)
    wcwt ~ Gamma(20.0,15.0/19)
    measerr ~ Gamma(10.0,20.0/9)
    staffweights ~ filldist(truncated(Normal(150,30),90.0,Inf),nstaff)
    patientweights ~ filldist(truncated(Normal(150,30),90.0,Inf),npatients)
    theta = LazyArray(@~ view(staffweights,staffid) .+ view(patientweights,patientid) .+ usewch .* wcwt)
    totweight ~ arraydist(LazyArray(@~ truncatenormal.(theta,measerr)))
end

But it gives a backtrace:

ERROR: MethodError: no method matching arraydist(::BroadcastVector{Any, typeof(truncatenormal), Tuple{BroadcastVector{Float64, typeof(+), Tuple{BroadcastVector{Float64, typeof(+), Tuple{SubArray{Float64, 1, Vector{Float64}, Tuple{Vector{Int64}}, false}, SubArray{Float64, 1, Vector{Float64}, Tuple{Vector{Int64}}, false}}}, BroadcastVector{Float64, typeof(*), Tuple{Vector{Int64}, Float64}}}}, Float64}})
Closest candidates are:
  arraydist(::AbstractVector{<:UnivariateDistribution}) at /home/dlakelan/.julia/packages/DistributionsAD/b93cZ/src/arraydist.jl:5
  arraydist(::AbstractMatrix{<:UnivariateDistribution}) at /home/dlakelan/.julia/packages/DistributionsAD/b93cZ/src/arraydist.jl:17
  arraydist(::AbstractVector{<:MultivariateDistribution}) at /home/dlakelan/.julia/packages/DistributionsAD/b93cZ/src/arraydist.jl:47
Stacktrace:
  [1] estweights3lazy(__model__::DynamicPPL.Model{typeof(estweights3lazy), (:nstaff, :staffid, :npatients, :patientid, :usewch, :totweight), (), (), Tuple{Int64, Vector{Int64}, Int64, Vector{Int64}, Vector{Int64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, __varinfo__::DynamicPPL.ThreadSafeVarInfo{DynamicPPL.UntypedVarInfo{DynamicPPL.Metadata{Dict{AbstractPPL.VarName, Int64}, Vector{Distribution}, Vector{AbstractPPL.VarName}, Vector{Real}, Vector{Set{DynamicPPL.Selector}}}, Float64}, Vector{Base.RefValue{Float64}}}, __context__::DynamicPPL.SamplingContext{DynamicPPL.SampleFromUniform, DynamicPPL.DefaultContext, Random._GLOBAL_RNG}, nstaff::Int64, staffid::Vector{Int64}, npatients::Int64, patientid::Vector{Int64}, usewch::Vector{Int64}, totweight::Vector{Float64})
    @ Main ./REPL[81]:7
  [2] macro expansion
    @ ~/.julia/packages/DynamicPPL/RcfQU/src/model.jl:465 [inlined]
  [3] _evaluate(model::DynamicPPL.Model{typeof(estweights3lazy), (:nstaff, :staffid, :npatients, :patientid, :usewch, :totweight), (), (), Tuple{Int64, Vector{Int64}, Int64, Vector{Int64}, Vector{Int64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, varinfo::DynamicPPL.ThreadSafeVarInfo{DynamicPPL.UntypedVarInfo{DynamicPPL.Metadata{Dict{AbstractPPL.VarName, Int64}, Vector{Distribution}, Vector{AbstractPPL.VarName}, Vector{Real}, Vector{Set{DynamicPPL.Selector}}}, Float64}, Vector{Base.RefValue{Float64}}}, context::DynamicPPL.SamplingContext{DynamicPPL.SampleFromUniform, DynamicPPL.DefaultContext, Random._GLOBAL_RNG})
    @ DynamicPPL ~/.julia/packages/DynamicPPL/RcfQU/src/model.jl:448
  [4] evaluate_threadsafe(model::DynamicPPL.Model{typeof(estweights3lazy), (:nstaff, :staffid, :npatients, :patientid, :usewch, :totweight), (), (), Tuple{Int64, Vector{Int64}, Int64, Vector{Int64}, Vector{Int64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, varinfo::DynamicPPL.UntypedVarInfo{DynamicPPL.Metadata{Dict{AbstractPPL.VarName, Int64}, Vector{Distribution}, Vector{AbstractPPL.VarName}, Vector{Real}, Vector{Set{DynamicPPL.Selector}}}, Float64}, context::DynamicPPL.SamplingContext{DynamicPPL.SampleFromUniform, DynamicPPL.DefaultContext, Random._GLOBAL_RNG})
...
82

Please post a full MWE for me to copy and paste.

83

Sure, it’s the same as the above MWE, but here it is in just copy-paste form:

using Pkg
Pkg.activate(".")
using Turing, DataFrames,DataFramesMeta,LazyArrays,Distributions,DistributionsAD
using LazyArrays

## every few hours a random staff member comes and gets a random
## patient to bring them outside to a garden through a door that has a
## scale. Sometimes using a wheelchair, sometimes not. knowing the
## total weight of the two people and the wheelchair plus some errors
## (from the scale measurements), infer the individual weights of all
## individuals and the weight of the wheelchair.

nstaff = 100
npat = 100
staffids = collect(1:nstaff)
patientids = collect(1:npat)
staffweights = rand(Normal(150,30),length(staffids))
patientweights = rand(Normal(150,30),length(staffids))
wheelchairwt = 15
nobs = 300

data = DataFrame(staff=rand(staffids,nobs),patient=rand(patientids,nobs))
data.usewch = rand(0:1,nobs)
data.totweights = [staffweights[data.staff[i]] + patientweights[data.patient[i]] for i in 1:nrow(data)] .+ data.usewch .* wheelchairwt .+ rand(Normal(0.0,20.0),nrow(data))


function truncatenormal(a,b)::UnivariateDistribution
    truncated(Normal(a,b),0.0,Inf)
end


@model function estweights3lazy(nstaff,staffid,npatients,patientid,usewch,totweight)
    wcwt ~ Gamma(20.0,15.0/19)
    measerr ~ Gamma(10.0,20.0/9)
    staffweights ~ filldist(truncated(Normal(150,30),90.0,Inf),nstaff)
    patientweights ~ filldist(truncated(Normal(150,30),90.0,Inf),npatients)
    theta = LazyArray(@~ view(staffweights,staffid) .+ view(patientweights,patientid) .+ usewch .* wcwt)
    totweight ~ arraydist(LazyArray(@~ truncatenormal.(theta,measerr)))
end

ch3l = sample(estweights3lazy(nstaff,data.staff,npat,data.patient,data.usewch,data.totweights),NUTS(500,.75),1000)
1 Like
84

I am not getting that error but I am getting numerical errors. I am using Turing 0.18.

85

as am I, let me just restart my Julia session and see if it persists.

Do you get it sampling at all?

86

also when I print the type of the BroadcastVector, I get

BroadcastVector{Truncated{Normal{ReverseDiff.TrackedReal{Float64, Float64, Nothing}}, Continuous, ReverseDiff.TrackedReal{Float64, Float64, Nothing}}, typeof(truncatenormal), Tuple{BroadcastVector{ReverseDiff.TrackedReal{Float64, Float64, Nothing}, typeof(+), Tuple{BroadcastVector{ReverseDiff.TrackedReal{Float64, Float64, Nothing}, typeof(+), Tuple{SubArray{ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}}, 1, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, Tuple{Vector{Int64}}, false}, SubArray{ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}}, 1, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, Tuple{Vector{Int64}}, false}}}, BroadcastVector{ReverseDiff.TrackedReal{Float64, Float64, Nothing}, typeof(*), Tuple{Vector{Int64}, ReverseDiff.TrackedReal{Float64, Float64, Nothing}}}}}, ReverseDiff.TrackedReal{Float64, Float64, Nothing}}}

which seems right to me.

Do you get it sampling at all?

No. I suspect the prior or model might be bad. Try a smaller dataset or initialising from a good point.

87

So since the model and data are entirely synthetic, the model is actually essentially perfect. All the real weights are normally distributed from the given priors and the true wheelchair weight is right at the peak density of the prior for the wcwt parameter. However obviously initialization might be the issue. I’ll see if I can figure out an init vector to avoid that… after I get my kids breakfast etc.

Thanks so much for your help!

1 Like
88

If it’s synthetic from a single set of parameter values and there is sufficient data and your model is identifiable, the posterior will asymptotically tend to a dirac-delta distribution. This will cause numerical issues when you have a large enough data set since you are trying to sample from a very concentrated distribution where the log joint almost everywhere is -Inf. A solution to this is to generate data from multiple sets of parameter values or use a smaller data set.

1 Like
89

Yes, or a very good initial point and a very small step size, good point.

However in this case there are 100 patients and 100 staff, and 300 observations so on average each person is only represented around 3 times, so I don’t think that is the issue, but I’m happy to try with 50 or 100 observations and see if it changes anything.

I let this model run for 1 hour and it was still at 0% and the backtrace when I interrupted looked like:

ERROR: InterruptException:
Stacktrace:
  [1] Array
    @ ./boot.jl:457 [inlined]
  [2] Array
    @ ./boot.jl:466 [inlined]
  [3] Array
    @ ./boot.jl:474 [inlined]
  [4] similar
    @ ./abstractarray.jl:829 [inlined]
  [5] similar
    @ ./abstractarray.jl:828 [inlined]
  [6] similar
    @ ./broadcast.jl:212 [inlined]
  [7] similar
    @ ./broadcast.jl:211 [inlined]
  [8] copy
    @ ./broadcast.jl:929 [inlined]
  [9] materialize
    @ ./broadcast.jl:904 [inlined]
 [10] logabsdetjac(b::Bijectors.TruncatedBijector{1, Float64, Float64}, x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:wcwt, :measerr, :staffweights, :patientweights), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:wcwt, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:wcwt, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:measerr, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:measerr, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:staffweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:staffweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:patientweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:patientweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(estweights3lazy), (:nstaff, :staffid, :npatients, :patientid, :usewch, :totweight), (), (), Tuple{Int64, Vector{Int64}, Int64, Vector{Int64}, Vector{Int64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 10}})
    @ Bijectors ~/.julia/packages/Bijectors/EELoe/src/bijectors/truncated.jl:125
 [11] logpdf_with_trans
    @ ~/.julia/packages/Bijectors/EELoe/src/Bijectors.jl:134 [inlined]
 [12] assume(rng::Random._GLOBAL_RNG, spl::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, dist::Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}, vn::AbstractPPL.VarName{:patientweights, Tuple{}}, vi::DynamicPPL.ThreadSafeVarInfo{DynamicPPL.TypedVarInfo{NamedTuple{(:wcwt, :measerr, :staffweights, :patientweights), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:wcwt, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:wcwt, Tuple{}}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:wcwt, :measerr, :staffweights, :patientweights), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:wcwt, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:wcwt, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:measerr, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:measerr, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:staffweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:staffweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:patientweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:patientweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(estweights3lazy), (:nstaff, :staffid, :npatients, :patientid, :usewch, :totweight), (), (), Tuple{Int64, Vector{Int64}, Int64, Vector{Int64}, Vector{Int64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 10}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:measerr, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:measerr, Tuple{}}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:wcwt, :measerr, :staffweights, :patientweights), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:wcwt, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:wcwt, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:measerr, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:measerr, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:staffweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:staffweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:patientweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:patientweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(estweights3lazy), (:nstaff, :staffid, :npatients, :patientid, :usewch, :totweight), (), (), Tuple{Int64, Vector{Int64}, Int64, Vector{Int64}, Vector{Int64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 10}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:staffweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:staffweights, Tuple{}}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:wcwt, :measerr, :staffweights, :patientweights), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:wcwt, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:wcwt, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:measerr, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:measerr, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:staffweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:staffweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:patientweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:patientweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(estweights3lazy), (:nstaff, :staffid, :npatients, :patientid, :usewch, :totweight), (), (), Tuple{Int64, Vector{Int64}, Int64, Vector{Int64}, Vector{Int64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 10}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:patientweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:patientweights, Tuple{}}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:wcwt, :measerr, :staffweights, :patientweights), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:wcwt, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:wcwt, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:measerr, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:measerr, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:staffweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:staffweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:patientweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:patientweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(estweights3lazy), (:nstaff, :staffid, :npatients, :patientid, :usewch, :totweight), (), (), Tuple{Int64, Vector{Int64}, Int64, Vector{Int64}, Vector{Int64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 10}}, Vector{Set{DynamicPPL.Selector}}}}}, ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:wcwt, :measerr, :staffweights, :patientweights), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:wcwt, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:wcwt, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:measerr, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:measerr, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:staffweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:staffweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:patientweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:patientweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(estweights3lazy), (:nstaff, :staffid, :npatients, :patientid, :usewch, :totweight), (), (), Tuple{Int64, Vector{Int64}, Int64, Vector{Int64}, Vector{Int64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 10}}, Vector{Base.RefValue{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:wcwt, :measerr, :staffweights, :patientweights), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:wcwt, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:wcwt, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:measerr, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:measerr, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:staffweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:staffweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:patientweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:patientweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(estweights3lazy), (:nstaff, :staffid, :npatients, :patientid, :usewch, :totweight), (), (), Tuple{Int64, Vector{Int64}, Int64, Vector{Int64}, Vector{Int64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 10}}}})
    @ Turing.Inference ~/.julia/packages/Turing/uMQmD/src/inference/hmc.jl:484
 [13] tilde_assume
    @ ~/.julia/packages/DynamicPPL/RcfQU/src/context_implementations.jl:63 [inlined]
 [14] tilde_assume
    @ ~/.julia/packages/DynamicPPL/RcfQU/src/context_implementations.jl:58 [inlined]
 [15] tilde_assume
    @ ~/.julia/packages/DynamicPPL/RcfQU/src/context_implementations.jl:43 [inlined]
 [16] tilde_assume!
    @ ~/.julia/packages/DynamicPPL/RcfQU/src/context_implementations.jl:140 [inlined]
 [17] estweights3lazy(__model__::DynamicPPL.Model{typeof(estweights3lazy), (:nstaff, :staffid, :npatients, :patientid, :usewch, :totweight), (), (), Tuple{Int64, Vector{Int64}, Int64, Vector{Int64}, Vector{Int64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, __varinfo__::DynamicPPL.ThreadSafeVarInfo{DynamicPPL.TypedVarInfo{NamedTuple{(:wcwt, :measerr, :staffweights, :patientweights), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:wcwt, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:wcwt, Tuple{}}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:wcwt, :measerr, :staffweights, :patientweights), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:wcwt, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:wcwt, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:measerr, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:measerr, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:staffweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:staffweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:patientweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:patientweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(estweights3lazy), (:nstaff, :staffid, :npatients, :patientid, :usewch, :totweight), (), (), Tuple{Int64, Vector{Int64}, Int64, Vector{Int64}, Vector{Int64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 10}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:measerr, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:measerr, Tuple{}}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:wcwt, :measerr, :staffweights, :patientweights), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:wcwt, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:wcwt, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:measerr, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:measerr, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:staffweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:staffweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:patientweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:patientweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(estweights3lazy), (:nstaff, :staffid, :npatients, :patientid, :usewch, :totweight), (), (), Tuple{Int64, Vector{Int64}, Int64, Vector{Int64}, Vector{Int64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 10}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:staffweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:staffweights, Tuple{}}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:wcwt, :measerr, :staffweights, :patientweights), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:wcwt, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:wcwt, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:measerr, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:measerr, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:staffweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:staffweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:patientweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:patientweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(estweights3lazy), (:nstaff, :staffid, :npatients, :patientid, :usewch, :totweight), (), (), Tuple{Int64, Vector{Int64}, Int64, Vector{Int64}, Vector{Int64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 10}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:patientweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:patientweights, Tuple{}}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:wcwt, :measerr, :staffweights, :patientweights), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:wcwt, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:wcwt, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:measerr, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:measerr, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:staffweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:staffweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:patientweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:patientweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(estweights3lazy), (:nstaff, :staffid, :npatients, :patientid, :usewch, :totweight), (), (), Tuple{Int64, Vector{Int64}, Int64, Vector{Int64}, Vector{Int64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 10}}, Vector{Set{DynamicPPL.Selector}}}}}, ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:wcwt, :measerr, :staffweights, :patientweights), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:wcwt, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:wcwt, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:measerr, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:measerr, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:staffweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:staffweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:patientweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:patientweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(estweights3lazy), (:nstaff, :staffid, :npatients, :patientid, :usewch, :totweight), (), (), Tuple{Int64, Vector{Int64}, Int64, Vector{Int64}, Vector{Int64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 10}}, Vector{Base.RefValue{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:wcwt, :measerr, :staffweights, :patientweights), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:wcwt, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:wcwt, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:measerr, Tuple{}}, Int64}, Vector{Gamma{Float64}}, Vector{AbstractPPL.VarName{:measerr, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:staffweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:staffweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:patientweights, Tuple{}}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64}, FillArrays.Fill{Truncated{Normal{Float64}, Continuous, Float64}, 1, Tuple{Base.OneTo{Int64}}}}}, Vector{AbstractPPL.VarName{:patientweights, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(estweights3lazy), (:nstaff, :staffid, :npatients, :patientid, :usewch, :totweight), (), (), Tuple{Int64, Vector{Int64}, Int64, Vector{Int64}, Vector{Int64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 10}}}}, __context__::DynamicPPL.SamplingContext{DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext, Random._GLOBAL_RNG}, nstaff::Int64, staffid::Vector{Int64}, npatients::Int64, patientid::Vector{Int64}, usewch::Vector{Int64}, totweight::Vector{Float64})

...

I tried 50 observations, and put the initial position at exactly the true values and the initial epsilon in NUTS at 1e-9,

as follows:

julia> ch3l = sample(estweights3lazy(nstaff,data.staff,npat,data.patient,data.usewch,data.totweights),NUTS(500,.75;init_ϵ=1e-9),1000;init_theta=vcat([15.0,20.0],staffweights,patientweights))

it still spits out a bunch of warnings like:

┌ Warning: The current proposal will be rejected due to numerical error(s).
│   isfinite.((θ, r, ℓπ, ℓκ)) = (true, true, false, true)
└ @ AdvancedHMC ~/.julia/packages/AdvancedHMC/HQHnm/src/hamiltonian.jl:47
┌ Warning: The current proposal will be rejected due to numerical error(s).
│   isfinite.((θ, r, ℓπ, ℓκ)) = (true, true, false, true)
└ @ AdvancedHMC ~/.julia/packages/AdvancedHMC/HQHnm/src/hamiltonian.jl:47
┌ Warning: The current proposal will be rejected due to numerical error(s).
│   isfinite.((θ, r, ℓπ, ℓκ)) = (true, true, false, true)
└ @ AdvancedHMC ~/.julia/packages/AdvancedHMC/HQHnm/src/hamiltonian.jl:47

and goes nowhere.

90

@cscherrer Would this be a good usecase for Soss.jl?

1 Like
91

Thanks @Akatz for the ping. There’s a lot here, is this the right place to start?

@model function fit300norms(norms,index)
    sd ~ Gamma(2.0,1.0)
    means ~ filldist(Uniform(0,50),300)
    for i in 1:300
        norms[index[i]] ~ Normal(means[index[i]],sd)
    end
end
92

I think we can rule out a model vs a computational issue by seeing if Stan will run this model. I suspect it will, and it will sample it nearly instantaneously. I personally suspect what’s going on is that ReverseDiff is failing to work properly. I’ll write up a Stan model and try it with Stan.jl, which will be a useful exercise for me anyway.

@csherrer, I think the MWE to start with is in this comment: Making Turing Fast with large numbers of parameters? - #51 by dlakelan

1 Like
93

This was kinda my feeling. I’m happy to make a PR, but I don’t currently understand Turing’s internals,
LazyArrays, or how they interact well enough to tell someone else in simple terms “if you have this kind of situation, use lazyarray inside arraydist (but not this other situation), and here’s why…”

My question was getting at this…but also whether it would be possible/desirable to have arraydist automatically convert its argument to a LazyArray.

I made a naive attempt at this, modifying the model from @Storopoli’s logistic regression example, and did get a sampling speedup relative to a simple arraydist (though not as much as with the bare LazyArray):

auto_lazy_arraydist(x) =  arraydist(BroadcastArray(xi -> xi, x))

@model logreg(X,  y; predictors=size(X, 2)) = begin
    #priors
    α ~ Normal(0, 2.5)
    β ~ filldist(TDist(3), predictors)

    #likelihood
    # Like this, samples in 50 s
    # y .~ BernoulliLogit.(α .+ X * β)

    # Like this, samples in 45 s
    # y ~ arraydist(BernoulliLogit.(α .+ X * β))

    # Like this, samples in 34 s
    y ~ auto_lazy_arraydist(BernoulliLogit.(α .+ X * β))

    # Like this, samples in 15 s
    # y ~ arraydist(LazyArray(@~ BernoulliLogit.(α .+ X * β)))
end;

X = Matrix(select(wells, Not(:switch)))
y = wells[:, :switch]
model = logreg(X, y);
chain = sample(model, NUTS(), MCMCThreads(), 2_000, 4)

Or maybe .~ could automatically parse to do a lazy broadcast or something? Ideally, I feel like users souldn’t need to learn about this external package in the first place, especially since some of its syntax (i.e., @~) is confusingly close to Turing’s ~, but with a completely different meaning.

Sorry if this comes across as greedy or demanding…but that is a risk when you are working with Julia programmers :wink:

1 Like
94

Perhaps some kind of macro defined in Turing:

@lazyarraydist(...broadcasted computation...)

which automatically makes a lazy array from a broadcasted array calculation

1 Like
95

Ugh, now I realize that I don’t have Stan installed, and this means I have to install Stan, which is a whole kettle of fish by itself … so I probably am not going to do this right away.

Since I got the arraydist with lazy arrays thing working ish, I am going to see if it works well with gamma distributions, and then compare to the non-working of the truncated normal. This is easier and faster to work on so it might produce some results more quickly.

Ok, so in fact the gamma version with lazy arrays and 300 observations samples quickly, a minute or so with a few numerical warnings.


@model function estweightslazygamma(nstaff,staffid,npatients,patientid,usewch,totweight)
    wcwt ~ Gamma(20.0,15.0/19)
    measerr ~ Gamma(10.0,20.0/9)
    staffweights ~ filldist(truncated(Normal(150,30),90.0,Inf),nstaff)
    patientweights ~ filldist(truncated(Normal(150,30),90.0,Inf),npatients)
    theta = LazyArray(@~ view(staffweights,staffid) .+ view(patientweights,patientid) .+ usewch .* wcwt)
    totweight ~ arraydist(LazyArray(@~ Gamma.(15, theta ./ 14)))
end

ch5 = sample(estweightslazygamma(nstaff,data.staff,npat,data.patient,data.usewch,data.totweights),NUTS(500,.75),1000)

So now I’m guessing that truncated distributions specifically somehow use if statements that prevent ReverseDiff from compiling the tape?

96

I didn’t use tape compilation when I used ReverseDiff above. You can generally rule out AD issues if you try 2 different reverse-mode packages and get similar behaviours.

97

Did you get a version using Truncated Normal to work? For me it always fails to do anything after tens of minutes to hours.

98

No but I didn’t try changing anything. When you get sampling issues, printing numbers inside the model can be helpful debugging info. Also try starting from just the prior, then add a few observations at a time and see where sampling stops working.

99

Well, debug printf let me realize I wasn’t turning on reverse diff, so inconsistent results depending on whether I was using a brand new Julia vs one where reverse diff had been enabled before…

debug printf also shows me that either it’s not really using my initialization or something weird is going on… I initialize it as:

ch3l = sample(estweights3lazy(nstaff,data.staff,npat,data.patient,data.usewch,data.totweights),NUTS(500,.75,init_ϵ=.002),1000;
              init_theta = vcat([15.0,20.0],staffweights,patientweights))

and then use this println:

    println("""Evaluating model... 
wcwt: $wcwt
measerr: $measerr
exstaffweights: $(staffweights[1:10])
expatweights: $(patientweights[1:10])
""")

inside the model

So initially wcwt and measerr should be 15.0 and 20.0, and all the estimated weights should be exactly correct… but the first print command says:

Evaluating model... 
wcwt: 0.6727279716358009
measerr: 0.19350227287063274
exstaffweights: [164.7901299487081, 217.52372933554884, 137.6205473016958, 175.4254441162211, 188.02642630005334, 175.8530682604024, 215.85351
585263436, 113.79632650916365, 109.91764961851418, 99.87309516255067]
expatweights: [168.51306657319324, 130.17891134846482, 148.83838569329237, 192.57929839739077, 137.26394389850617, 149.17636016275125, 174.854
40785493637, 109.5090910975809, 151.17248643226876, 205.4490091028178]

so measerr goes to nearly zero and of course this means serious numerical issues. Why is the initial vector not being used?

100

could be a bug, please open an issue

1 Like