-
Notifications
You must be signed in to change notification settings - Fork 218
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Problems with deterministic distribution #1335
Comments
Isn't the 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 Ultimatively, I don't think any of these ideas will yield a perfect solution since it seems you need rather a solution to TuringLang/DynamicPPL.jl#94 (the general problem in your example is a duplicate of this issue IMO). Letting HMC (or any other sampler) deal with So a maybe better solution for now would be to just perform these computations afterwards based on your samples of |
BTW regardless of your problem with @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
|
Thanks! I was just trying to replicate I tried your |
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 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 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 |
For anyone checking out this issue now, the |
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 👍 |
@torfjelde Maybe add this to the Turing docs? |
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 👍 |
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 |
Sounds good to me too if that is possible. |
I was only referring to the first part of your comment, actually. Couldn't the |
I would have to check, if the string macros already allow this then even better. |
Weell, last time I checked, So I think we need more than just defining a more generic |
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 Maybe I miss something, but to me it seems |
That's true; my point is more that because the separation of the packages has been quite aggressive, the
Yeah, exactly. That's what I meant with:
in the above 👍 |
I thought more, we could change |
Ah, yeah we could just iterate through the chains and concatenate 👍 |
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 |
Could we add the following to DynamicPPL (or something similar that works, it is untested 😛)? 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 |
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
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,)
... |
This PR adds `generated_quantities` as discussed in TuringLang/Turing.jl#1335 + adds a fix for #167.
I am also trying to return derived parameters from a Turing.jl Is there an easy way to return a multivariate deterministic variable? (aka a vector) |
Related TuringLang/DynamicPPL.jl#94 |
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:
Typical results look like this:
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:
The text was updated successfully, but these errors were encountered: