Turing.jl - Binomial ArgumentError zero(p) <= p <= one(p) is not satisfied

mangulo · April 15, 2020, 12:17am

Hi All,

This is my first post in Julia discourse! I have been trying to do parameter inference in a classical Susceptible-Infected-Recovered epidemic model using Turing.
I am finding a strange error: the program correctly runs sometimes, but in other cases it returns the error:

ERROR: ArgumentError: Binomial: the condition zero(p) <= p <= one(p) is not satisfied.

Any suggestion about what is going wrong here? Code to reproduce the error is pasted below. Thanks a lot!

##  Generate data

n = 8.556e6 # total population
S0 = 8.556e6
I0 = 1
R0 = 0

β = 0.6
γ = 1/5

number_iterations = 70

S = Array{Int64}(undef, number_iterations)
I = Array{Int64}(undef, number_iterations)
R = Array{Int64}(undef, number_iterations)

S[1] = S0; I[1] = I0; R[1] = R0

@inbounds for t = 2 :number_iterations

    # Transitions between states
    si = Binomial(S[t-1],  β* I[t-1] / n) |> rand
    ir = Binomial(I[t-1], γ) |> rand

    # Equations
    S[t] = S[t-1] - si
    I[t] = I[t-1] + si - ir
    R[t] = R[t-1] + ir

end

S_true = deepcopy(S)
I_true = deepcopy(I)
R_true = deepcopy(R)

plot(I_true)

@model SIR_model_POMP(I, R) = begin
        n = 8.556e6 # total population

        β  ~ truncated(Normal(0.7, 0.2), 0, 1)
        γ  ~ truncated(Normal(1/5,0.01), 1/6, 1/4)

        #  Number of observations.
        number_iterations = length(S)

        S[1] = 8.556e6 # initial susceptible population

        for t = 2:number_iterations
                # Transitions
                si ~ Binomial(S[t-1], β* I[t-1] / n)
                ir ~ Binomial(I[t-1], γ)

                # Equations
                S[t] = S[t-1] - si
                I[t] = I[t-1] + si - ir
                R[t] = R[t-1] + ir
        end
end


g = Gibbs(HMC(1e-4, 7, :β, :γ))
model = SIR_model_POMP(I_true, R_true)
tchain = mapreduce(c -> sample(model, g, 4000), chainscat, 1:4);

trappmartin · April 15, 2020, 9:50am

@Kai_Xu do you have an idea why this happens?

mohamed82008 · April 15, 2020, 2:41pm

I think your p parameter is either more than 1 or less than 0 by a small amount due to numerical error. Try scaling it between 0 and 1 before constructing the Binomial distribution constructor.

Edit: it might also be out of range because of a problem in your model so print out the value in the model and see if it’s out of range by a small or a large amount. A large amount is a sign that something else is going on.

mangulo · April 15, 2020, 3:55pm

Thanks a lot for your suggestion! I modified the code following your suggestion, but I am still getting the same error:

ERROR: ArgumentError: Binomial: the condition zero(p) <= p <= one(p) is not satisfied.

I also added the array temparray to track the value of the p argument given to the Binomial function. It looks like this before the code crashes:

julia> temparray[1:141796]
141796-element Array{Any,1}:
 9.714446093935619e-8 
 1.9428892187871238e-7
 2.9143338281806855e-7
 3.8857784375742476e-7
 4.857223046967809e-7 
 5.828667656361371e-7 
 6.800112265754933e-7 
 7.771556875148495e-7 
 8.743001484542057e-7 
 9.714446093935619e-7 
 1.068589070332918e-6 
 1.1657335312722742e-6
 ⋮                    
 0.0                  
 0.0                  
 0.0                  
 0.0                  
 0.0                  
 0.0                  
 0.0                  
 0.0                  
 0.0                  
 0.0                  
 0.0                  
 0.0

And at the step 141797, the value of p given to the Binomial function looks like this:

julia> temparray[141797]
Dual{ForwardDiff.Tag{Turing.Core.var"#f#7"{DynamicPPL.VarInfo{NamedTuple{(:β, :γ, :si, :ir),Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:β},Int64},Array{Truncated{Normal{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:β},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:γ},Int64},Array{Truncated{Normal{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:γ},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:si},Int64},Array{Binomial{Float64},1},Array{DynamicPPL.VarName{:si},1},Array{Int64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:ir},Int64},Array{Binomial{Float64},1},Array{DynamicPPL.VarName{:ir},1},Array{Int64,1},Array{Set{DynamicPPL.Selector},1}}}},Float64},DynamicPPL.Model{var"##inner_function#4367#452",NamedTuple{(:I, :R),Tuple{Array{Int64,1},Array{Int64,1}}},DynamicPPL.ModelGen{(:I, :R),var"###SIR_model_POMP#4393",NamedTuple{(),Tuple{}}},Val{()}},DynamicPPL.Sampler{HMC{Turing.Core.ForwardDiffAD{40},(:β, :γ),AdvancedHMC.UnitEuclideanMetric},Turing.Inference.HMCState{DynamicPPL.VarInfo{NamedTuple{(:β, :γ, :si, :ir),Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:β},Int64},Array{Truncated{Normal{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:β},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:γ},Int64},Array{Truncated{Normal{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:γ},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:si},Int64},Array{Binomial{Float64},1},Array{DynamicPPL.VarName{:si},1},Array{Int64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:ir},Int64},Array{Binomial{Float64},1},Array{DynamicPPL.VarName{:ir},1},Array{Int64,1},Array{Set{DynamicPPL.Selector},1}}}},Float64},AdvancedHMC.StaticTrajectory{AdvancedHMC.EndPointTS,AdvancedHMC.Leapfrog{Float64}},AdvancedHMC.Adaptation.NoAdaptation,AdvancedHMC.PhasePoint{Array{Float64,1},AdvancedHMC.DualValue{Float64,Array{Float64,1}}}}}},Float64}}(8.747967470967713e-8,2.2003229232244862e-8,0.0)

And when the code crashes it looks like this

julia> temparray[end]
Dual{ForwardDiff.Tag{Turing.Core.var"#f#7"{DynamicPPL.VarInfo{NamedTuple{(:β, :γ, :si, :ir),Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:β},Int64},Array{Truncated{Normal{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:β},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:γ},Int64},Array{Truncated{Normal{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:γ},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:si},Int64},Array{Binomial{Float64},1},Array{DynamicPPL.VarName{:si},1},Array{Int64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:ir},Int64},Array{Binomial{Float64},1},Array{DynamicPPL.VarName{:ir},1},Array{Int64,1},Array{Set{DynamicPPL.Selector},1}}}},Float64},DynamicPPL.Model{var"##inner_function#4410#466",NamedTuple{(:I, :R),Tuple{Array{Int64,1},Array{Int64,1}}},DynamicPPL.ModelGen{(:I, :R),var"###SIR_model_POMP#4436",NamedTuple{(),Tuple{}}},Val{()}},DynamicPPL.Sampler{HMC{Turing.Core.ForwardDiffAD{40},(:β, :γ),AdvancedHMC.UnitEuclideanMetric},Turing.Inference.HMCState{DynamicPPL.VarInfo{NamedTuple{(:β, :γ, :si, :ir),Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:β},Int64},Array{Truncated{Normal{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:β},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:γ},Int64},Array{Truncated{Normal{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:γ},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:si},Int64},Array{Binomial{Float64},1},Array{DynamicPPL.VarName{:si},1},Array{Int64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:ir},Int64},Array{Binomial{Float64},1},Array{DynamicPPL.VarName{:ir},1},Array{Int64,1},Array{Set{DynamicPPL.Selector},1}}}},Float64},AdvancedHMC.StaticTrajectory{AdvancedHMC.EndPointTS,AdvancedHMC.Leapfrog{Float64}},AdvancedHMC.Adaptation.NoAdaptation,AdvancedHMC.PhasePoint{Array{Float64,1},AdvancedHMC.DualValue{Float64,Array{Float64,1}}}}}},Float64}}(NaN,NaN,NaN)

Any further suggestion of what is going wrong? Thanks!

The revised code is pasted below:

##  Generate data

n = 8.556e6 # total population
S0 = 8.556e6
I0 = 1
R0 = 0

β = 0.6
γ = 1/5

number_iterations = 70

S = Array{Int64}(undef, number_iterations)
I = Array{Int64}(undef, number_iterations)
R = Array{Int64}(undef, number_iterations)

S[1] = S0; I[1] = I0; R[1] = R0

@inbounds for t = 2 :number_iterations

    # Transitions between states
    si = Binomial(S[t-1],  β* I[t-1] / n) |> rand
    ir = Binomial(I[t-1], γ) |> rand

    # Equations
    S[t] = S[t-1] - si
    I[t] = I[t-1] + si - ir
    R[t] = R[t-1] + ir

end

S_true = deepcopy(S)
I_true = deepcopy(I)
R_true = deepcopy(R)

#plot(I_true)



## Inference

global temparray = []

clip(p) = minimum([maximum([0, p]), 1])
ReLU(p) = maximum([0, p])

@model SIR_model_POMP(I, R) = begin
        n = 8.556e6 # total population

        β  ~ truncated(Normal(0.7, 0.2), 0, 1)
        γ  ~ truncated(Normal(1/5,0.01), 1/6, 1/4)

        #  Number of observations.
        number_iterations = length(S)

        S0 = 8.556e6
        I0 = 1
        R0 = 0

        S[1] = 8.556e6 # initial susceptible population
        I[1] = I0
        R[1] = R0

        for t = 2:number_iterations
                # Transitions

                temp = clip(β*I[t-1]/n)
                push!(temparray, temp)

                si ~ Binomial(S[t-1], temp)
                ir ~ Binomial(I[t-1], γ)

                # Equations
                S[t] = ReLU(S[t-1] - si)
                I[t] = ReLU(I[t-1] + si - ir)
                R[t] = ReLU(R[t-1] + ir)
        end
end


g = Gibbs(HMC(1e-4, 7, :β, :γ))
model = SIR_model_POMP(I_true, R_true)
tchain = mapreduce(c -> sample(model, g, 10), chainscat, 1:10);

mohamed82008 · April 15, 2020, 4:19pm

I the way you are using Turing is wrong. The lines:

si ~ Binomial(S[t-1], temp)
ir ~ Binomial(I[t-1], γ)

tell Turing that si and ir are random variables or parameters to be sampled. These random variables are sampled more than once in the loop which is not allowed in Turing and then they are not assigned any sampler in the Gibbs sampler constructor. The reason for the error is that you get NaN for β for some reason, most likely because of some interplay between the above mistakes leading to unpredictable results.

mohamed82008 · April 15, 2020, 4:29pm

Arguably we should be throwing an error with your example though explaining how to fix it, so please open an issue

mangulo · April 15, 2020, 4:57pm

Thanks, @mohamed82008! I understand your point.
What would be the correct way to implement the above SIR epidemic model in Turing?
Namely, the discrete time equations are:

S_{t+1} = S_{t} - s_t, \quad I_{t+1} = I_t + s_t - r_t, \quad R_{t+1} = R_t + r_t,

where S_t, I_t and R_t denote the number of Susceptible, Infected and Recovered individuals at time t, respectively.

Above, s_t and r_t are random variables representing the new infected and the new recovered individuals at time t, respectively. These random variables satisfy s_t \sim \mbox{Binomial}(S_t, \beta I_t /n) and r_t \sim \mbox{Binomial}(I_t, \gamma). The parameters of the SIR model are \beta and \gamma.

In Julia, the SIR model would read like this:

##  Generate data

n = 8.556e6 # total population
S0 = 8.556e6
I0 = 1
R0 = 0

β = 0.6
γ = 1/5

number_iterations = 70

S = Array{Int64}(undef, number_iterations)
I = Array{Int64}(undef, number_iterations)
R = Array{Int64}(undef, number_iterations)

S[1] = S0; I[1] = I0; R[1] = R0

@inbounds for t = 2 :number_iterations

    # Transitions between states
    si = Binomial(S[t-1],  β* I[t-1] / n) |> rand
    ir = Binomial(I[t-1], γ) |> rand

    # Equations
    S[t] = S[t-1] - si
    I[t] = I[t-1] + si - ir
    R[t] = R[t-1] + ir

end

Thanks!

mangulo · April 15, 2020, 9:48pm

Not sure how to do this. Do you have a link explaining the process of opening an issue? Thanks!

mohamed82008 · April 15, 2020, 11:15pm

The first thing you should define is what the observations are. ir and si seem to be reasonable observations, i.e. then number of people recovered and infected every day. So you should make those vectors and pass them in the model arguments:

@model SIR_model_POMP(ir, si)
    ...
end

Then inside the loop, do ir[t] ~ ... and si[t] ~ .... This will treat ir and si as observations which I think makes sense in this context.

On opening issues: Creating an issue - GitHub Docs.

mangulo · April 17, 2020, 2:26am

Thanks for your suggestion. Certainly, the number of recovered ir and infected si individuals can be observed in the simple SIR epidemiological model. However, for the applicability of Turing to other more general epidemiological models, it is more natural to keep the total the total number of infected I and recovered R individuals as observables.

In that sense, it would be better if one could use a Turing model of the following form:

@model SIR_model_POMP(I, R) = begin
        n = 8.556e6 # total population

        β  ~ truncated(Normal(0.7, 0.2), 0, 1)
        γ  ~ truncated(Normal(1/5,0.01), 1/6, 1/4)

        #  Number of observations.
        number_iterations = length(I)

        S0 = 8.556e6
        I0 = 1
        R0 = 0

        S, I, R = SIR_model(S0, I0, R0, n, β, γ, number_iterations)
end

In this way, one would avoid the sampling problem introduced by the variables si and ir, and the subsequent problems you mention.

Above, the SIR_model function returns the vectors S, I, and R with the total number of Susceptible, Infected and Recovered individuals given its input parameters. Namely, it reads like this:

function SIR_model(S0, I0, R0, n, β, γ, number_iterations)

        S = Array{Int64}(undef, number_iterations)
        I = Array{Int64}(undef, number_iterations)
        R = Array{Int64}(undef, number_iterations)

        S[1] = S0; I[1] = I0; R[1] = R0

        for t = 2 :number_iterations

            # Transitions between states
            si = Binomial(S[t-1],  β* I[t-1] / n) |> rand
            ir = Binomial(I[t-1], γ) |> rand

            # Equations
            S[t] = S[t-1] - si
            I[t] = I[t-1] + si - ir
            R[t] = R[t-1] + ir

        end
        return S, I, R
end

Finally, one would like to call a sampler to infer the parameters β and γ:

g = Gibbs(HMC(1e-4, 7, :β, :γ))
model = SIR_model_POMP(I_true, R_true)
tchain = mapreduce(c -> sample(model, g, 10), chainscat, 1:10);

With the above code I get the error:

ERROR: MethodError: no method matching Distributions.BinomialGeomSampler(::Int64, ::ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#7"{DynamicPPL.VarInfo{NamedTuple{(:β, :γ),Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:β},Int64},Array{Truncated{Normal{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:β},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:γ},Int64},Array{Truncated{Normal{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:γ},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}}}},Float64},DynamicPPL.Model{var"##inner_function#503#10",NamedTuple{(:I, :R),Tuple{Array{Int64,1},Array{Int64,1}}},DynamicPPL.ModelGen{(:I, :R),var"###SIR_model_POMP#517",NamedTuple{(),Tuple{}}},Val{()}},DynamicPPL.Sampler{HMC{Turing.Core.ForwardDiffAD{40},(:β, :γ),AdvancedHMC.UnitEuclideanMetric},Turing.Inference.SamplerState{DynamicPPL.VarInfo{NamedTuple{(:β, :γ),Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:β},Int64},Array{Truncated{Normal{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:β},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:γ},Int64},Array{Truncated{Normal{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:γ},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}}}},Float64}}}},Float64},Float64,2})

Seems to me that Turing does not know how to differentiate the SIR_model function.
Is there some way to make a code like this work in Turing?
Thanks!

I paste below the code to reproduce the error.

##  Generate data
function SIR_model(S0, I0, R0, n, β, γ, number_iterations)

        S = Array{Int64}(undef, number_iterations)
        I = Array{Int64}(undef, number_iterations)
        R = Array{Int64}(undef, number_iterations)

        S[1] = S0; I[1] = I0; R[1] = R0

        for t = 2 :number_iterations

            # Transitions between states
            si = Binomial(S[t-1],  β* I[t-1] / n) |> rand
            ir = Binomial(I[t-1], γ) |> rand

            # Equations
            S[t] = S[t-1] - si
            I[t] = I[t-1] + si - ir
            R[t] = R[t-1] + ir

        end
        return S, I, R
end


n = 8.556e6 # total population
S0 = 8.556e6
I0 = 1
R0 = 0

β = 0.6
γ = 1/5

number_iterations = 70

S_true, I_true, R_true = SIR_model(S0, I0, R0, n, β, γ, number_iterations)

## Inference

@model SIR_model_POMP(I, R) = begin
        n = 8.556e6 # total population

        β  ~ truncated(Normal(0.7, 0.2), 0, 1)
        γ  ~ truncated(Normal(1/5,0.01), 1/6, 1/4)

        #  Number of observations.
        number_iterations = length(I)

        S0 = 8.556e6
        I0 = 1
        R0 = 0

        S, I, R = SIR_model(S0, I0, R0, n, β, γ, number_iterations)
end


g = Gibbs(HMC(1e-4, 7, :β, :γ))
model = SIR_model_POMP(I_true, R_true)
tchain = mapreduce(c -> sample(model, g, 10), chainscat, 1:10);