cscherrer / Soss.jl

Probabilistic programming via source rewriting
https://cscherrer.github.io/Soss.jl/stable/
MIT License
413 stars 30 forks source link

Can't get particles from a model with a Poisson likelihood #99

Open sethaxen opened 4 years ago

sethaxen commented 4 years ago

Calling particles on a model with a Poisson distributed variable fails if its rate parameter is drawn from another distribution. eg:

julia> using Soss, Distributions

julia> m = @model begin
           λ ~ Normal(10000, 1)
           y ~ Poisson(λ)
       end;

julia> particles(m())
ERROR: UndefVarError: maybe_particles not defined
Stacktrace:
 [1] <(::Particles{Float64,1000}, ::Int64) at /Users/saxen/.julia/packages/MonteCarloMeasurements/21tbg/src/register_primitive.jl:45
 [2] sampler at /Users/saxen/.julia/packages/Distributions/KjaXI/src/univariate/discrete/poisson.jl:147 [inlined]
 [3] rand at /Users/saxen/.julia/packages/Distributions/KjaXI/src/univariates.jl:158 [inlined]
 [4] rand at /Users/saxen/.julia/packages/Distributions/KjaXI/src/genericrand.jl:26 [inlined]
 [5] #Particles#18(::Bool, ::Bool, ::Type{Particles}, ::Random._GLOBAL_RNG, ::Int64, ::Poisson{Particles{Float64,1000}}) at /Users/saxen/.julia/packages/MonteCarloMeasurements/21tbg/src/types.jl:63
 [6] Particles at /Users/saxen/.julia/packages/MonteCarloMeasurements/21tbg/src/types.jl:60 [inlined]
 [7] #Particles#19 at /Users/saxen/.julia/packages/MonteCarloMeasurements/21tbg/src/types.jl:68 [inlined]
 [8] Particles at /Users/saxen/.julia/packages/MonteCarloMeasurements/21tbg/src/types.jl:68 [inlined]
 [9] parts at /Users/saxen/projects/Soss.jl/src/particles.jl:41 [inlined]
 [10] macro expansion at /Users/saxen/.julia/packages/GeneralizedGenerated/EBwdX/src/closure_conv.jl:121 [inlined]
 [11] _particles(::Type{TypeEncoding(Main)}, ::Model{NamedTuple{(),T} where T<:Tuple,TypeEncoding(begin
    λ ~ Normal(10000, 1)
    y ~ Poisson(λ)
end),TypeEncoding(Main)}, ::NamedTuple{(),Tuple{}}, ::Val{1000}) at /Users/saxen/.julia/packages/GeneralizedGenerated/EBwdX/src/closure_conv.jl:121
 [12] particles at /Users/saxen/projects/Soss.jl/src/particles.jl:66 [inlined] (repeats 2 times)
 [13] top-level scope at REPL[3]:1
cscherrer commented 4 years ago

If we unwind what Soss is doing, we get to

julia> λ = Particles(Normal(10000,1))
Part500(10000.0 ± 1.0)

julia> pois = Poisson(λ)
Poisson{Particles{Float64,500}}(
λ: 10000.0 ± 1.0
)

julia> Particles(500, pois)
ERROR: UndefVarError: maybe_particles not defined

Originally all of the parts and particles weirdness in Soss was to hack around some temporary limitations of MonteCarloMeasurements, while avoiding type piracy. @baggepinnen is there a more direct way to do this sort of thing now?

baggepinnen commented 4 years ago

I can not reproduce it here, which version of MonteCarloMeasurements are you running?

julia> λ = Particles(Normal(10000,1))
Part500(10000.0 ± 1.0)

julia> pois = Poisson(λ)
Poisson{Particles{Float64,500}}(
λ: 10000.0 ± 1.0
)

julia> Particles(500, pois)
ERROR: TypeError: non-boolean (Particles{Float64,500}) used in boolean context
Stacktrace:
 [1] sampler at /home/fredrikb/.julia/packages/Distributions/KjaXI/src/univariate/discrete/poisson.jl:147 [inlined]
 [2] rand at /home/fredrikb/.julia/packages/Distributions/KjaXI/src/univariates.jl:158 [inlined]
 [3] rand at /home/fredrikb/.julia/packages/Distributions/KjaXI/src/genericrand.jl:24 [inlined]
 [4] #Particles#110(::Bool, ::Bool, ::Type{Particles}, ::Int64, ::Poisson{Particles{Float64,500}}) at /home/fredrikb/.julia/packages/MonteCarloMeasurements/hYwGr/src/particles.jl:238
 [5] Particles(::Int64, ::Poisson{Particles{Float64,500}}) at /home/fredrikb/.julia/packages/MonteCarloMeasurements/hYwGr/src/particles.jl:235
 [6] top-level scope at none:0

I tried on Soss v0.9.0 which installs MCM v0.5.8

EDIT: I could reproduce the problem on master. Should be solved by the commit referenced below.

baggepinnen commented 4 years ago

I think I have found and solved the problem in https://github.com/baggepinnen/MonteCarloMeasurements.jl/commit/84d5fb2c2930cdf362fcbe524bf84642d4fa5119

baggepinnen commented 4 years ago

A new version is tagged, I hope it resolves this issue.

sethaxen commented 4 years ago

Thanks @baggepinnen, the fix resolved that issue! I'm now seeing the error you posted:

julia> using Soss, Distributions
[ Info: Precompiling Soss [8ce77f84-9b61-11e8-39ff-d17a774bf41c]

julia> m = @model begin
           λ ~ Normal(10000, 1)
           y ~ Poisson(λ)
       end;

julia> particles(m())
ERROR: TypeError: non-boolean (Particles{Bool,1000}) used in boolean context
Stacktrace:
 [1] sampler at /Users/saxen/.julia/packages/Distributions/KjaXI/src/univariate/discrete/poisson.jl:147 [inlined]
 [2] rand at /Users/saxen/.julia/packages/Distributions/KjaXI/src/univariates.jl:158 [inlined]
 [3] rand at /Users/saxen/.julia/packages/Distributions/KjaXI/src/genericrand.jl:26 [inlined]
 [4] #Particles#18(::Bool, ::Bool, ::Type{Particles}, ::Random._GLOBAL_RNG, ::Int64, ::Poisson{Particles{Float64,1000}}) at /Users/saxen/.julia/packages/MonteCarloMeasurements/6Onbj/src/types.jl:63
 [5] Particles at /Users/saxen/.julia/packages/MonteCarloMeasurements/6Onbj/src/types.jl:60 [inlined]
 [6] #Particles#19 at /Users/saxen/.julia/packages/MonteCarloMeasurements/6Onbj/src/types.jl:68 [inlined]
 [7] Particles at /Users/saxen/.julia/packages/MonteCarloMeasurements/6Onbj/src/types.jl:68 [inlined]
 [8] parts at /Users/saxen/projects/Soss.jl/src/particles.jl:41 [inlined]
 [9] macro expansion at /Users/saxen/.julia/packages/GeneralizedGenerated/EBwdX/src/closure_conv.jl:121 [inlined]
 [10] _particles(::Type{TypeEncoding(Main)}, ::Model{NamedTuple{(),T} where T<:Tuple,TypeEncoding(begin
    λ ~ Normal(10000, 1)
    y ~ Poisson(λ)
end),TypeEncoding(Main)}, ::NamedTuple{(),Tuple{}}, ::Val{1000}) at /Users/saxen/.julia/packages/GeneralizedGenerated/EBwdX/src/closure_conv.jl:121
 [11] particles at /Users/saxen/projects/Soss.jl/src/particles.jl:66 [inlined] (repeats 2 times)
 [12] top-level scope at REPL[3]:1

This arises when the rate parameter of the Poisson distribution is a Particles and this boolean is checked: https://github.com/JuliaStats/Distributions.jl/blob/master/src/univariate/discrete/poisson.jl#L147

@cscherrer, @baggepinnen what's the recommended way to proceed here? This is the issue with having a Distribution with Particles parameters instead of a Particles{Distribution} right?

baggepinnen commented 4 years ago

Yeah that question is a bit harder to answer, see the full discussion here https://github.com/baggepinnen/MonteCarloMeasurements.jl/issues/22

I have a draft solution that operates on the representation Distribution{Particles}, but think I have some compelling arguments in favor of the representation Particles{Distribution} instead. When you want the former is essentially when you create or display the distribution, but for computations, the latter representation is more convenient.