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

**URL:** <https://discourse.julialang.org/t/turing-jl-binomial-argumenterror-zero-p-p-one-p-is-not-satisfied/37601>\
**Category:** Probabilistic Programming\
**Tags:** turing\
**Created:** [April 15, 2020, 12:17am UTC](https://discourse.julialang.org/t/turing-jl-binomial-argumenterror-zero-p-p-one-p-is-not-satisfied/37601 "2020-04-15T00:17:34Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![mangulo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mangulo/32/14076_2.png) [@mangulo](https://discourse.julialang.org/u/mangulo)\
**Post date:** [April 15, 2020, 12:17am UTC](https://discourse.julialang.org/t/turing-jl-binomial-argumenterror-zero-p-p-one-p-is-not-satisfied/37601/1 "2020-04-15T00:17:34Z")

</div>

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!

```julia
## 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);

```

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [April 15, 2020, 9:50am UTC](https://discourse.julialang.org/t/turing-jl-binomial-argumenterror-zero-p-p-one-p-is-not-satisfied/37601/2 "2020-04-15T09:50:04Z")

</div>

@Kai_Xu do you have an idea why this happens?

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [April 15, 2020, 2:41pm UTC](https://discourse.julialang.org/t/turing-jl-binomial-argumenterror-zero-p-p-one-p-is-not-satisfied/37601/3 "2020-04-15T14:41:53Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mangulo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mangulo/32/14076_2.png) [@mangulo](https://discourse.julialang.org/u/mangulo)\
**Post date:** [April 15, 2020, 3:55pm UTC](https://discourse.julialang.org/t/turing-jl-binomial-argumenterror-zero-p-p-one-p-is-not-satisfied/37601/4 "2020-04-15T15:55:30Z")

</div>

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
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:

```julia
## 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);

```

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [April 15, 2020, 4:19pm UTC](https://discourse.julialang.org/t/turing-jl-binomial-argumenterror-zero-p-p-one-p-is-not-satisfied/37601/5 "2020-04-15T16:19:18Z")

</div>

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [April 15, 2020, 4:29pm UTC](https://discourse.julialang.org/t/turing-jl-binomial-argumenterror-zero-p-p-one-p-is-not-satisfied/37601/6 "2020-04-15T16:29:12Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mangulo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mangulo/32/14076_2.png) [@mangulo](https://discourse.julialang.org/u/mangulo)\
**Post date:** [April 15, 2020, 4:57pm UTC](https://discourse.julialang.org/t/turing-jl-binomial-argumenterror-zero-p-p-one-p-is-not-satisfied/37601/7 "2020-04-15T16:57:59Z")

</div>

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:

```julia
## 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!

---

<div class="post-metadata">

**Author:** ![mangulo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mangulo/32/14076_2.png) [@mangulo](https://discourse.julialang.org/u/mangulo)\
**Post date:** [April 15, 2020, 9:48pm UTC](https://discourse.julialang.org/t/turing-jl-binomial-argumenterror-zero-p-p-one-p-is-not-satisfied/37601/8 "2020-04-15T21:48:27Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [April 15, 2020, 11:15pm UTC](https://discourse.julialang.org/t/turing-jl-binomial-argumenterror-zero-p-p-one-p-is-not-satisfied/37601/9 "2020-04-15T23:15:44Z")

</div>

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:

```julia
@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](https://help.github.com/en/github/managing-your-work-on-github/creating-an-issue).

---

<div class="post-metadata">

**Author:** ![mangulo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mangulo/32/14076_2.png) [@mangulo](https://discourse.julialang.org/u/mangulo)\
**Post date:** [April 17, 2020, 2:26am UTC](https://discourse.julialang.org/t/turing-jl-binomial-argumenterror-zero-p-p-one-p-is-not-satisfied/37601/10 "2020-04-17T02:26:33Z")

</div>

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:

```julia
@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:

```julia
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 `γ`:

```julia
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.

```julia
## 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);

```