Turing.jl: Is Binomial supported/working?

[goedman]

Hi, I'm trying to get below Turing program to run:

using Turing, MCMCChain

data = (n = 9, k = 6)
@model globe_toss(n, k) = begin
  #prior
  theta ~ Beta(1, 1)
  #model
  k ~ Binomial(n, theta)
end

chn2 = sample(globe_toss(data.n, data.k), HMC(2000, 0.75, 5))
describe(chn2)

which produces below error:

julia> include("/Users/rob/.julia/dev/StatisticalRethinking/chapters/02/turing_globe_toss.jl")
[ Info: Assume - theta is a parameter
[ Info: Observe - k is an observation
ERROR: LoadError: MethodError: no method matching Float64(::Flux.Tracker.TrackedReal{Float64})
Closest candidates are:
Float64(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:185
Float64(::T<:Number) where T<:Number at boot.jl:725
Float64(::Int8) at float.jl:60
...

This is on Julia 1.0.2 and using Turing v0.5.1. These are my first steps in using Turing, so it is likely I'm doing something wrong. Thanks.

[yebai]

This should be already supported. We have a similar model here:
https://github.com/TuringLang/TuringExamples/blob/master/stan-models/bernoulli.model.jl

Ps. seems like we run into a Flux.Tracker issue here.

[cpfiffer]

@yebai I've noticed this come up a couple times while I've been playing with some of the models. It seems like the parameters are being converted into tracked values which don't shoehorn well into the binomlogpdf function. I suspect this might need more triage because I've seen this come up a good deal.

The way I've patched this is by calling Flux.Tracker.data() on the parameters. It's not a good long-term solution but it certainly will get things working. Here's the modified version of @goedman's example:

using Turing, MCMCChain, Flux.Tracker

data = (n = 9, k = 6)
@model globe_toss(n, k) = begin
  #prior
  theta ~ Beta(1, 1)
  #model
  k ~ Binomial(n, Tracker.data(theta))
end

chn2 = sample(globe_toss(data.n, data.k), HMC(2000, 0.75, 5))
describe(chn2)

[yebai]

k ~ Binomial(n, Tracker.data(theta))

@cpfiffer
Unfortunately, this hack will break the correctness of AD and gradient-based samplers. Could you give it a try to AD through binomlogpdf using FLux/ForwardDiff outside Turing? I suspect Distributions.binomlogpdf is not correctly implemented to support Dual number types from AD packages such as Flux.

[xukai92]

I wonder maybe it's a Flux.jl issue or some change in Distributions.jl. There are pdf functions of generic support defined in StatsFuns.jl. But they are not correctly trigger in the current version.

[goedman]

Thank you @cpfiffer, that works well! The results look correct for this toy example (compared with CmdStan and e.g. grid approximation). For my purpose that is better than falling back to Bernoulli.

@yebai does your statement imply the above example can also break? Can I use it until the Flux issue is resolved?

Some background: I'm trying to 'translate' the StatisticalRethinking book to Julia and I'm looking into using the Turing PPL as a replacement for map() in McElreath's R package rethinking:

globe2.qa <- map(
    alist(
        w ~ dbinom(9, p) ,  # binomial likelihood
        p ~ dbeta(1, 1)       # beta(1, 1) prior
    ) ,
    data=list(w=6) )

# display summary of quadratic approximation
precis( globe2.qa )

In addition to map() (which will be the @model name(...) begin .... end construct) I also want to create a mapping to CmdStan equivalent to McElreath's map2stan().

The other option would be to attempt to follow a brms route that uses modeling like MixedModels.jl,
but my preference is to follow the book as closely as possible because that is what makes this book so special (hence my remark about Bernoulli above).

[willtebbutt]

@goedman could you try running your original example on the master branch of Turing? Hopefully #601 has resolved this.

[goedman]

Thank you! Fixed.

[cpfiffer]

@willtebbutt Do you know if this is related? This model seems to be failing in either :forward_diff or :reverse_diff mode:

using Turing
using StatsFuns: logistic

Turing.setadbackend(:reverse_diff)

@model m10_4(y, actors, x₁, x₂) = begin
   # Number of unique actors in the data set
   N_actor = length(unique(actors))

   # Total number of data points for all actors
   # actor = tzeros(length(actors))

   # Array which’ll contain priors for each unique actor
   α = TArray{Any}(undef, N_actor)
   println(typeof(α))

   # For each actor [1...7] set a prior
   for i ∈ 1:length(α)
       α[i] ~ Normal(0,10)
   end

   βp ~ Normal(0, 10)
   βpC ~ Normal(0, 10)

   for i ∈ 1:length(y)
       p = logistic(α[actors[i]] + (βp + βpC * x₁[i]) * x₂[i])

       y[i] ~ Binomial(1, p)
   end
end

z = sample(m10_4([1.0,2.0,3.0], [1,1,2], [1,1,1], [2,2,2]), HMC(300, 0.01, 5))

Output:

ERROR: LoadError: MethodError: no method matching Float64(::Flux.Tracker.TrackedReal{Float64})
Closest candidates are:
  Float64(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:185
  Float64(::T<:Number) where T<:Number at boot.jl:725
  Float64(::Int8) at float.jl:60
  ...
Stacktrace:
 [1] convert(::Type{Float64}, ::Flux.Tracker.TrackedReal{Float64}) at ./number.jl:7
 [2] cconvert(::Type, ::Flux.Tracker.TrackedReal{Float64}) at ./essentials.jl:344
 [3] binomlogpdf(::Int64, ::Flux.Tracker.TrackedReal{Float64}, ::Int64) at /home/cameron/.julia/packages/StatsFuns/0W2sM/src/rmath.jl:58
 [4] logpdf(::Binomial{Flux.Tracker.TrackedReal{Float64}}, ::Int64) at /home/cameron/.julia/packages/Distributions/WHjOk/src/univariates.jl:526
 [5] logpdf at /home/cameron/.julia/packages/Distributions/WHjOk/src/univariates.jl:325 [inlined]
 [6] macro expansion at ./logging.jl:308 [inlined]
...

[yebai]

@cpfiffer This should be working after upgrading Bijectors. See https://github.com/TuringLang/Bijectors.jl/pull/11