This example has been auto-generated from the examples/ folder at GitHub repository.

Coin toss model (Beta-Bernoulli)

# Activate local environment, see `Project.toml`
import Pkg; Pkg.activate(".."); Pkg.instantiate();

In this example, we are going to perform an exact inference for a coin toss model that can be represented as:

\[\begin{aligned} p(\theta) &= \mathrm{Beta}(\theta|a, b),\\ p(y_i|\theta) &= \mathrm{Ber}(y_i|\theta),\\ \end{aligned}\]

where $y_i \in \{0, 1\}$ is a binary observation induced by Bernoulli likelihood while $\theta$ is a Beta prior distribution on the parameter of Bernoulli. We are interested in inferring the posterior distribution of $\theta$.

We start with importing all needed packages:

using RxInfer, Random

Let's generate some synthetic dataset with IID observations from Bernoulli distribution, that represents our coin tosses. We also assume that our coin is biased:

rng = MersenneTwister(42)
n = 500
θ_real = 0.75
distribution = Bernoulli(θ_real)

dataset = float.(rand(rng, Bernoulli(θ_real), n));

# GraphPPL.jl export `@model` macro for model specification
# It accepts a regular Julia function and builds an FFG under the hood
@model function coin_model(n)

    # `datavar` creates data 'inputs' in our model
    # We will pass data later on to these inputs
    # In this example we create a sequence of inputs that accepts Float64
    y = datavar(Float64, n)

    # We endow θ parameter of our model with some prior
    θ ~ Beta(4.0, 8.0)
    # or, in this particular case, the `Uniform(0.0, 1.0)` prior also works:
    # θ ~ Uniform(0.0, 1.0)

    # We assume that outcome of each coin flip is governed by the Bernoulli distribution
    for i in 1:n
        y[i] ~ Bernoulli(θ)
    end

end

result = inference(
    model = coin_model(length(dataset)), 
    data  = (y = dataset, )
)

Inference results:
  Posteriors       | available for (θ)

θestimated = result.posteriors[:θ]

Beta{Float64}(α=365.0, β=147.0)

using Plots

rθ = range(0, 1, length = 1000)

p = plot(title = "Inference results")

plot!(rθ, (x) -> pdf(Beta(2.0, 7.0), x), fillalpha=0.3, fillrange = 0, label="P(θ)", c=1,)
plot!(rθ, (x) -> pdf(θestimated, x), fillalpha=0.3, fillrange = 0, label="P(θ|y)", c=3)
vline!([θ_real], label="Real θ")