Turing.jl
Turing.jl copied to clipboard
TuringLang
Problems with deterministic distribution
I'm trying to use deterministic distribution in my model as outlined here: https://discourse.julialang.org/t/is-there-a-turing-alternative-to-pm-deterministic-from-pymc3/38667/2 and so far my code fails:
Distributions.rand(rng::AbstractRNG, d::Determin) = d.val
Distributions.logpdf(d::Determin, x::T) where T<:Real = zero(x)
Bijectors.bijector(d::Determin) = Identity{0}() # <= `0` indicates 0-dimensional, i.e. univariate
@model model1(a, ::Type{T} = Float64) where {T} = begin
ma ~ Normal(0.0, 0.5)
for i in 1:length(a)
a[i] ~ Normal(ma, 0.05)
end
b ~ Determin(-exp(-2*ma))
return ma, b
end
function test_model1()
# sampler = HMC(0.001, 10)
# sampler = HMCDA(1000, 0.75, 0.05)
sampler = NUTS(0.65)
c1 = sample(model1(a1), sampler, 2000)
println(c1)
end
Typical results look like this:
julia> test_model1()
┌ Info: Found initial step size
└ ϵ = 0.025
Object of type Chains, with data of type 1000×14×1 Array{Float64,3}
Iterations = 1:1000
Thinning interval = 1
Chains = 1
Samples per chain = 1000
internals = acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, is_accept, log_density, lp, max_hamiltonian_energy_error, n_steps, nom_step_size, numerical_error, step_size, tree_depth
parameters = b, ma
2-element Array{ChainDataFrame,1}
Summary Statistics
parameters mean std naive_se mcse ess r_hat
────────── ───────────────────────── ──────────────────────── ─────────────────────── ──────────────────────── ──────── ──────
b 13816403112283807744.0000 7133491167269512192.0000 225580797572648256.0000 2172001283591027712.0000 4.8969 1.3982
ma 0.0690 0.0179 0.0006 0.0007 830.2122 0.9998
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
────────── ───────────────────────── ───────────────────────── ───────────────────────── ───────────────────────── ─────────────────────────
b -5026632495806943232.0000 12090906044107700224.0000 15653120955249981440.0000 18317459011682740224.0000 22252367610115112960.0000
ma 0.0349 0.0565 0.0694 0.0811 0.1047
HMC and HMCDA samplers also give similarly bad results bud SMC actually works. Did I make something stupid here or is this a bug? A similar model in PyMC3 works with the NUTS sampler:
a1 = [0.05, 0.03, 0.02, 0.03, 0.04, 0.09, 0.13, 0.16]
with pm.Model() as model3:
mu_a = pm.Normal('mu_a', 0.0, 0.5)
aim = pm.Normal('aim', mu=mu_a, sigma=0.05, observed=a1)
b = pm.Deterministic('b', -np.exp(-2*mu_a))
tr = pm.sample(2000, tune=700)
pm.traceplot(tr)
display(pm.summary(tr))
| mean | sd | hdi_3% | hdi_97% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat
-- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | --
0.069 | 0.017 | 0.037 | 0.101 | 0.000 | 0.0 | 1860.0 | 1860.0 | 1851.0 | 2777.0 | 1.0
-0.872 | 0.030 | -0.923 | -0.811 | 0.001 | 0.0 | 1860.0 | 1860.0 | 1851.0 | 2777.0 | 1.0
Isn't the logpdf implementation incorrect and should rather be something like
Distributions.logpdf(d::Determin, x) = x == d.T ? zero(float(x)) : oftype(float(x), -Inf)
? I haven't performed any debugging yet but I assume that the incorrect logpdf will lead HMC to return samples that do not adhere to the restriction that you want to impose by the Determin distribution. I'm wondering as well if maybe a DiscreteNonParametric distribution with support only at a single point would be more appropriate. Alternatively, you could try to use a Normal distribution with zero variance and check what happens instead of implementing your own distributions and bijectors.
Ultimatively, I don't think any of these ideas will yield a perfect solution since it seems you need rather a solution to https://github.com/TuringLang/DynamicPPL.jl/issues/94 (the general problem in your example is a duplicate of this issue IMO). Letting HMC (or any other sampler) deal with b doesn't seem to be correct or helpful in general, you would rather want to just perform the computation AND track it.
So a maybe better solution for now would be to just perform these computations afterwards based on your samples of mu_a, and not during inference.
BTW regardless of your problem with Determin you can simplify (and probably speed up) your model by implementing it as
@model function model1(a)
ma ~ Normal(0.0, 0.5)
a ~ filldist(Normal(ma, 0.05), length(a)) # alternatively: a .~ Normal(ma, 0.05)
b ~ Determin(-exp(-2*ma))
return ma, b
end
T is never used and hence not needed, and filldist should be faster than the loop.
Thanks! I was just trying to replicate Deterministic from PyMC3. It's actually just a small part of my full model where T is actually used but I didn't know about filldist.
I tried your logpdf and it doesn't work either. For now I'll just keep using PyMC3 and I'll port my model when Turing is more mature.
I quickly ran your model with
b ~ Normal(-exp(-2*ma), 0)
instead of
b ~ Determin(-exp(-2*ma))
and as I assumed basically all samples of the HMC sampler are rejected (which makes inference impossible) due to the fact that the Hamiltonian dynamics don't respect the deterministic constraint. HMC doesn't seem to be suitable for sampling b (which ideally you wouldn't do at all...) - PG seems to work fine.
So I guess for now performing these computations on the resulting chain instead of inside the model might be the easiest option, if one wants to use Turing + HMC.
And to make what @devmotion suggested easier, you can potentially make use of the generated_quantities method from here: https://github.com/cambridge-mlg/Covid19/blob/4d967d58eeada3f3399116fc815cb17fd4163c87/src/utils.jl#L79.
With that, your code would end up something like:
@model model1(a, ::Type{T} = Float64) where {T} = begin
ma ~ Normal(0.0, 0.5)
for i in 1:length(a)
a[i] ~ Normal(ma, 0.05)
end
b = -exp(-2*ma)
return (b = b, ) # or just `b` if you don't want a named tuple
end
chain = sample(model, ...)
results = generated_quantities(model, chain) # <= results in an array of named tuples
You can also use the arrtup2tuparr method in the file I linked above to convert results into a slightly more convenient format.
For anyone checking out this issue now, the generated_quantities method is dependent on https://github.com/TuringLang/DynamicPPL.jl/pull/147 to be merged to work with the most recent versions of Turing.
Update:
function generated_quantities(m::Turing.Model, c::MCMCChains.Chains)
# if `c` is multiple chains we pool them into a single chain
chain = length(chains(c)) == 1 ? c : MCMCChains.pool_chain(c)
varinfo = Turing.DynamicPPL.VarInfo(m)
return map(1:length(chain)) do i
Turing.DynamicPPL.setval!(varinfo, chain, i, 1)
m(varinfo)
end
end
should now work nicely for all cases in the most recent release :+1:
I'd prefer to just add it to one of the packages, but it's a bit unclear which package it should go in. Once we figure that out, we could add it to the docs :+1:
Hi, I think it would be great if we could support this functionality in Turing. I suppose Turing.jl is the best place for it. Moreover, I would suggest adding functionality to obtain the log prob values of each datum at each iteration after inference. I am not aware that we have this atm. This is useful for model selection, see: https://github.com/TuringLang/MCMCChains.jl#deviance-information-criterion-dic, and would be a great feature especially if it can easily be used for other models too.
Maybe I'm missing something here, but doesn't generated_quantities above provide this functionality? So you are wondering where to put it? IMO it could also be put into DynamicPPL (which does not depend on MCMCChains though) or MCMCChains (which does not depend on DynamicPPL though) since it is not Turing.jl-specific - it doesn't depend on any objects or functionality in Turing, but only on stuff that belongs to MCMCChains and DynamicPPL. Maybe it would be possible to define it for AbstractChains and put it into DynamicPPL (DynamicPPL already depends on AbstractMCMC)? It seems the main thing would be to replace pool_chain with something more generic.
Sounds good to me too if that is possible. Regarding log prob, I’m aware that this should be possible I was merely asking to add an additional interface function for users that does that.
Regarding log prob, I’m aware that this should be possible I was merely asking to add an additional interface function for users that does that.
I was only referring to the first part of your comment, actually. Couldn't the prob/logprob macros be used to compute these log probabilities from the chains? Are they lacking some functionality?
I would have to check, if the string macros already allow this then even better. @mohamed82008
Maybe it would be possible to define it for AbstractChains and put it into DynamicPPL (DynamicPPL already depends on AbstractMCMC)? It seems the main thing would be to replace pool_chain with something more generic.
Weell, last time I checked, setval! and others implicitly depends on MCMCChains.jl by making certain assumptions about the indexing behaviour of AbstractMCMC.AbstractChains (which is not enforced/informed in AbstractMCMC.jl) :confused:
So I think we need more than just defining a more generic pool_chains method, tbh. Plus, the pool_chains could just be dropped and we add an optional chain_idx argument and inform the user that if they have several chains and want to run generated_quantities, use MCMCChains.pool_chains. Again, there's an implicit dependency, so it's not nice :confused:
Weell, last time I checked,
setval!and others implicitly depends onMCMCChains.jlby making certain assumptions about the indexing behaviour ofAbstractMCMC.AbstractChains(which is not enforced/informed inAbstractMCMC.jl) confused
I think that's fine for now, and actually I'm a bit reluctant with enforcing such things in AbstractMCMC in a top-down approach - I think it is better to test and use methods in downstream packages in this case and move them upstream if they have been shown to be useful, and similarly only make things part of the general interface if needed. It seems, the current implementation in https://github.com/TuringLang/DynamicPPL.jl/blob/4432591d2e149858cb59a6993a11038acaf71dc5/src/varinfo.jl#L1150 just depends on indexing with three indices and keys being defined. That's not too many assumptions, so I'm not too worried about it right now.
Maybe I miss something, but to me it seems pool_chains is not needed at all since the implementation of setval! in DynamicPPL already allows to specify the chains index as third argument - so you would just iterate through the chain indices as well instead of pooling the chains first and then always using the chain index 1, as you currently do.
That's not too many assumptions, so I'm not too worried about it right now.
That's true; my point is more that because the separation of the packages has been quite aggressive, the AbstractChains is really just a proxy for MCMCChains.Chains so that we can avoid direct dependence on AbstractChains.jl (I guess due to compilation times?). I don't have a suggestion for a better approach though; just noting that it's, uhmm, a bit unfortunate :sweat_smile:
Maybe I miss something, but to me it seems pool_chains is not needed at all since the implementation of setval! in DynamicPPL already allows to specify the chains index as third argument - so you would just iterate through the chain indices as well instead of pooling the chains first and then always using the chain index 1, as you currently do.
Yeah, exactly. That's what I meant with:
Plus, the pool_chains could just be dropped and we add an optional chain_idx argument and inform the user that if they have several chains and want to run generated_quantities, use MCMCChains.pool_chains.
in the above :+1:
I thought more, we could change generated_quantities such that this is done automatically (so users don't have to care about it).
Yes, or even return results as a matrix instead (corresponding to the structure of Chains).
Btw, seems like this is now the working impl:
function generated_quantities(m::Turing.Model, c::MCMCChains.Chains)
# if `c` is multiple chains we pool them into a single chain
chain = length(chains(c)) == 1 ? c : MCMCChains.pool_chain(c)
varinfo = Turing.DynamicPPL.VarInfo(m)
return map(1:length(chain)) do i
empty!(varinfo) # <= NEW
Turing.DynamicPPL.setval!(varinfo, chain, i, 1)
m(varinfo)
end
end
@devmotion you know why we now need to clear the varinfo to recompute the return-value? Before the return-value would be recomputed every no matter, I believe. I'm not suggesting that we change this, I'm just curious which change made it so.
Could we add the following to DynamicPPL (or something similar that works, it is untested :stuck_out_tongue:)?
function generated_quantities(model::Model, chain::AbstractChains)
varinfo = VarInfo(model)
iters = Iterators.product(1:size(chain, 1), 1:size(chain, 3))
return map(iters) do (sample_idx, chain_idx)
setval!(varinfo, chain, sample_idx, chain_idx)
model(varinfo)
end
end
We use basically the same logic for the logprob macro already. I'm not sure why empty! would be needed here - actually I am very surprised that it is (do you have a MWE?) since it is not used in the logprob macro implementation.
Hmm, now I'm super-confused. The original bug was that the values would be the same for all samples (even though the parameters were changing). Now I'm getting the following on the exact same problem with exactly the same package versions:
julia> using Turing
julia> function generated_quantities(m::Turing.Model, c::MCMCChains.Chains)
# if `c` is multiple chains we pool them into a single chain
chain = length(chains(c)) == 1 ? c : MCMCChains.pool_chain(c)
varinfo = DynamicPPL.VarInfo(m)
return map(1:length(chain)) do i
# empty!(varinfo)
DynamicPPL.setval!(varinfo, chain[i])
m(varinfo)
end
end
generated_quantities (generic function with 1 method)
julia> # Stuff
C = [0.089398 0.0267295 0.007992; 0.0 0.000164962 4.31057e-5; 0.0170703 0.00223028 0.000291394; 0.0 0.0 0.00431887]
4×3 Array{Float64,2}:
0.089398 0.0267295 0.007992
0.0 0.000164962 4.31057e-5
0.0170703 0.00223028 0.000291394
0.0 0.0 0.00431887
julia> @model dyadRel(C, ::Type{T}=Vector{Float64}) where {T} = begin
nCrosses = size(C)[2]
σ ~ InverseGamma(1, 10)
T2 = typeof(σ)
δ = T(undef, nCrosses)
for c in 1:nCrosses
δ[c] ~ truncated(Cauchy(0, σ), T2(0), T2(Inf))
end
Δ = δ / sum(δ)
lp = sum(log.(C * Δ))
Turing.acclogp!(_varinfo, lp)
θ = (sum(Δ[1])/2 + Δ[2]/4)
r = 2 * θ
return (r = r, )
end
┌ Warning: you are using the internal variable `_varinfo`
└ @ DynamicPPL ~/.julia/packages/DynamicPPL/moP7G/src/compiler.jl:169
dyadRel (generic function with 2 methods)
julia> model = dyadRel(C)
DynamicPPL.Model{var"#17#18",(:C, :T),(:T,),(),Tuple{Array{Float64,2},Type{Array{Float64,1}}},Tuple{Type{Array{Float64,1}}}}(var"#17#18"(), NamedTuple{(:C, :T),Tuple{Array{Float64,2},Type{Array{Float64,1}}}}(([0.089398 0.0267295 0.007992; 0.0 0.000164962 4.31057e-5; 0.0170703 0.00223028 0.000291394; 0.0 0.0 0.00431887], Array{Float64,1})), NamedTuple{(:T,),Tuple{Type{Array{Float64,1}}}}((Array{Float64,1},)))
julia> fit = sample(model, NUTS(0.65), 100);
┌ Info: Found initial step size
└ ϵ = 1.6
julia> res = generated_quantities(model, fit)
ERROR: BoundsError: attempt to access 1×16×1 Array{Float64,3} at index [16:16, 1:16, 1:1]
Stacktrace:
[1] throw_boundserror(::Array{Float64,3}, ::Tuple{UnitRange{Int64},Base.Slice{Base.OneTo{Int64}},Base.Slice{Base.OneTo{Int64}}}) at ./abstractarray.jl:541
[2] checkbounds at ./abstractarray.jl:506 [inlined]
[3] _getindex at ./multidimensional.jl:742 [inlined]
[4] getindex at ./abstractarray.jl:1060 [inlined]
[5] getindex at /home/tor/.julia/packages/AxisArrays/IFpjG/src/indexing.jl:104 [inlined]
[6] getindex(::Chains{Float64,AxisArrays.AxisArray{Float64,3,Array{Float64,3},Tuple{AxisArrays.Axis{:iter,StepRange{Int64,Int64}},AxisArrays.Axis{:var,Array{Symbol,1}},AxisArrays.Axis{:chain,UnitRange{Int64}}}},Missing,NamedTuple{(:parameters, :internals),Tuple{Array{Symbol,1},Array{Symbol,1}}},NamedTuple{(),Tuple{}}}, ::Int64, ::Function, ::Function) at /home/tor/.julia/packages/MCMCChains/8AK5L/src/chains.jl:113
[7] getindex at /home/tor/.julia/packages/MCMCChains/8AK5L/src/chains.jl:104 [inlined]
[8] #76 at /home/tor/.julia/packages/DynamicPPL/moP7G/src/varinfo.jl:1187 [inlined]
[9] mapreduce_first at ./reduce.jl:392 [inlined]
[10] _mapreduce(::DynamicPPL.var"#76#79"{Chains{Float64,AxisArrays.AxisArray{Float64,3,Array{Float64,3},Tuple{AxisArrays.Axis{:iter,StepRange{Int64,Int64}},AxisArrays.Axis{:var,Array{Symbol,1}},AxisArrays.Axis{:chain,UnitRange{Int64}}}},Missing,NamedTuple{(:parameters, :internals),Tuple{Array{Symbol,1},Array{Symbol,1}}},NamedTuple{(),Tuple{}}}}, ::typeof(vcat), ::IndexLinear, ::Array{Int64,1}) at ./reduce.jl:403
[11] _mapreduce_dim at ./reducedim.jl:318 [inlined]
[12] #mapreduce#620 at ./reducedim.jl:310 [inlined]
[13] mapreduce at ./reducedim.jl:310 [inlined]
[14] _setval_kernel!(::DynamicPPL.VarInfo{NamedTuple{(:σ, :δ),Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:σ,Tuple{}},Int64},Array{InverseGamma{Float64},1},Array{DynamicPPL.VarName{:σ,Tuple{}},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:δ,Tuple{Tuple{Int64}}},Int64},Array{Truncated{Cauchy{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:δ,Tuple{Tuple{Int64}}},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}}}},Float64}, ::DynamicPPL.VarName{:σ,Tuple{}}, ::Chains{Float64,AxisArrays.AxisArray{Float64,3,Array{Float64,3},Tuple{AxisArrays.Axis{:iter,StepRange{Int64,Int64}},AxisArrays.Axis{:var,Array{Symbol,1}},AxisArrays.Axis{:chain,UnitRange{Int64}}}},Missing,NamedTuple{(:parameters, :internals),Tuple{Array{Symbol,1},Array{Symbol,1}}},NamedTuple{(),Tuple{}}}, ::Array{Symbol,1}) at /home/tor/.julia/packages/DynamicPPL/moP7G/src/varinfo.jl:1186
[15] macro expansion at /home/tor/.julia/packages/DynamicPPL/moP7G/src/varinfo.jl:1169 [inlined]
[16] _typed_setval!(::DynamicPPL.VarInfo{NamedTuple{(:σ, :δ),Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:σ,Tuple{}},Int64},Array{InverseGamma{Float64},1},Array{DynamicPPL.VarName{:σ,Tuple{}},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:δ,Tuple{Tuple{Int64}}},Int64},Array{Truncated{Cauchy{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:δ,Tuple{Tuple{Int64}}},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}}}},Float64}, ::NamedTuple{(:σ, :δ),Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:σ,Tuple{}},Int64},Array{InverseGamma{Float64},1},Array{DynamicPPL.VarName{:σ,Tuple{}},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:δ,Tuple{Tuple{Int64}}},Int64},Array{Truncated{Cauchy{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:δ,Tuple{Tuple{Int64}}},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}}}}, ::Chains{Float64,AxisArrays.AxisArray{Float64,3,Array{Float64,3},Tuple{AxisArrays.Axis{:iter,StepRange{Int64,Int64}},AxisArrays.Axis{:var,Array{Symbol,1}},AxisArrays.Axis{:chain,UnitRange{Int64}}}},Missing,NamedTuple{(:parameters, :internals),Tuple{Array{Symbol,1},Array{Symbol,1}}},NamedTuple{(),Tuple{}}}, ::Array{Symbol,1}) at /home/tor/.julia/packages/DynamicPPL/moP7G/src/varinfo.jl:1160
[17] _setval! at /home/tor/.julia/packages/DynamicPPL/moP7G/src/varinfo.jl:1159 [inlined]
[18] setval! at /home/tor/.julia/packages/DynamicPPL/moP7G/src/varinfo.jl:1148 [inlined]
[19] #15 at ./REPL[24]:7 [inlined]
[20] iterate at ./generator.jl:47 [inlined]
[21] _collect(::UnitRange{Int64}, ::Base.Generator{UnitRange{Int64},var"#15#16"{DynamicPPL.Model{var"#17#18",(:C, :T),(:T,),(),Tuple{Array{Float64,2},Type{Array{Float64,1}}},Tuple{Type{Array{Float64,1}}}},Chains{Float64,AxisArrays.AxisArray{Float64,3,Array{Float64,3},Tuple{AxisArrays.Axis{:iter,StepRange{Int64,Int64}},AxisArrays.Axis{:var,Array{Symbol,1}},AxisArrays.Axis{:chain,UnitRange{Int64}}}},Missing,NamedTuple{(:parameters, :internals),Tuple{Array{Symbol,1},Array{Symbol,1}}},NamedTuple{(),Tuple{}}},DynamicPPL.VarInfo{NamedTuple{(:σ, :δ),Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:σ,Tuple{}},Int64},Array{InverseGamma{Float64},1},Array{DynamicPPL.VarName{:σ,Tuple{}},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:δ,Tuple{Tuple{Int64}}},Int64},Array{Truncated{Cauchy{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:δ,Tuple{Tuple{Int64}}},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}}}},Float64}}}, ::Base.EltypeUnknown, ::Base.HasShape{1}) at ./array.jl:699
[22] collect_similar(::UnitRange{Int64}, ::Base.Generator{UnitRange{Int64},var"#15#16"{DynamicPPL.Model{var"#17#18",(:C, :T),(:T,),(),Tuple{Array{Float64,2},Type{Array{Float64,1}}},Tuple{Type{Array{Float64,1}}}},Chains{Float64,AxisArrays.AxisArray{Float64,3,Array{Float64,3},Tuple{AxisArrays.Axis{:iter,StepRange{Int64,Int64}},AxisArrays.Axis{:var,Array{Symbol,1}},AxisArrays.Axis{:chain,UnitRange{Int64}}}},Missing,NamedTuple{(:parameters, :internals),Tuple{Array{Symbol,1},Array{Symbol,1}}},NamedTuple{(),Tuple{}}},DynamicPPL.VarInfo{NamedTuple{(:σ, :δ),Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:σ,Tuple{}},Int64},Array{InverseGamma{Float64},1},Array{DynamicPPL.VarName{:σ,Tuple{}},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:δ,Tuple{Tuple{Int64}}},Int64},Array{Truncated{Cauchy{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:δ,Tuple{Tuple{Int64}}},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}}}},Float64}}}) at ./array.jl:628
[23] map(::Function, ::UnitRange{Int64}) at ./abstractarray.jl:2162
[24] generated_quantities(::DynamicPPL.Model{var"#17#18",(:C, :T),(:T,),(),Tuple{Array{Float64,2},Type{Array{Float64,1}}},Tuple{Type{Array{Float64,1}}}}, ::Chains{Float64,AxisArrays.AxisArray{Float64,3,Array{Float64,3},Tuple{AxisArrays.Axis{:iter,StepRange{Int64,Int64}},AxisArrays.Axis{:var,Array{Symbol,1}},AxisArrays.Axis{:chain,UnitRange{Int64}}}},Missing,NamedTuple{(:parameters, :internals),Tuple{Array{Symbol,1},Array{Symbol,1}}},NamedTuple{(),Tuple{}}}) at ./REPL[24]:5
[25] top-level scope at REPL[29]:1
but it works if I add uncomment the empty! statement.
We use basically the same logic for the logprob macro already.
This is indeed where I originally got the idea from:)
Using your impl @devmotion, I get the error again (I dunno..):
julia> function generated_quantities(model::DynamicPPL.Model, chain::MCMCChains.Chains)
varinfo = DynamicPPL.VarInfo(model)
iters = Iterators.product(1:size(chain, 1), 1:size(chain, 3))
return map(iters) do (sample_idx, chain_idx)
DynamicPPL.setval!(varinfo, chain, sample_idx, chain_idx)
model(varinfo)
end
end
generated_quantities (generic function with 2 methods)
julia> res = generated_quantities(model, fit)
50×1 Array{NamedTuple{(:r,),Tuple{Float64}},2}:
(r = 0.3297491812718108,)
(r = 0.3297491812718108,)
(r = 0.3297491812718108,)
(r = 0.3297491812718108,)
(r = 0.3297491812718108,)
...
but
julia> function generated_quantities(model::DynamicPPL.Model, chain::MCMCChains.Chains)
varinfo = DynamicPPL.VarInfo(model)
iters = Iterators.product(1:size(chain, 1), 1:size(chain, 3))
return map(iters) do (sample_idx, chain_idx)
empty!(varinfo)
DynamicPPL.setval!(varinfo, chain, sample_idx, chain_idx)
model(varinfo)
end
end
generated_quantities (generic function with 2 methods)
julia> res = generated_quantities(model, fit)
50×1 Array{NamedTuple{(:r,),Tuple{Float64}},2}:
(r = 0.4816375116174122,)
(r = 0.7905255445789598,)
(r = 0.6521361515488626,)
(r = 0.49979176287674854,)
(r = 0.7775244794392783,)
(r = 0.7989344174804943,)
(r = 0.46010138685964236,)
...
I am also trying to return derived parameters from a Turing.jl @model but finding it extremely difficult.
Is there an easy way to return a multivariate deterministic variable? (aka a vector)
