Quick start

The code examples below introduce different modeling frameworks, or–if you know what you’re looking for–get you going right away with a specific model type.

Logistic regression

Have a look at this most simple logistic regression!

R

library(brms)

n <- 30
x <- rnorm(n, mean = 0, sd = 1)
p <- ifelse(x > 1, 0.2, 0.8)
y <- sapply(p, function(p) sample(c(0, 1), 1, prob = c(1 - p, p)))

data <- data.frame(x = x, y = y)

fit <- brm(y ~ x, family = bernoulli(link = "logit"), data = data)
fit

Python

import numpy as np
import pymc as pm
import arviz as az

N = 30
x = np.random.normal(0, 1, N)
p = np.where(x > 1, 0.2, 0.8)
y = np.random.binomial(1, p)

with pm.Model() as model:
    alpha = pm.Normal("alpha", mu=0, sigma=1)
    beta = pm.Normal("beta", mu=0, sigma=1)

    mu = alpha + beta * x

    p = pm.math.invlogit(mu)
    y_obs = pm.Bernoulli("y_obs", p=p, observed=y)

    idata = pm.sample(1000, tune=1000, chains=4, return_inferencedata=True)

summary = az.summary(idata)

Julia

using Turing, StatsFuns

N = 30
x = rand(Normal(0, 1), N)
p = [i > 1 ? 0.2 : 0.8 for i in x]
y = rand.(Bernoulli.(p))

@model function demo(x, y)
    α ~ Normal(0, 1)
    β ~ Normal(0, 1)

    for i in eachindex(y)
        p = logistic(α + β * x[i])
        y[i] ~ Bernoulli(p)
    end
end

model = demo(x, y)
chain = sample(model, NUTS(), MCMCThreads(), 1000, 4)

Multilevel models Work In Progress

Hierarchical modeling comes natural in the Bayesian framework. Simply stack distributions!

Let’s first have a look at the simplest case: a random intercept-only model.

R

#

Python

#

Julia

using Turing, LinearAlgebra

@model function random_intercept(group, y)
    n_groups = length(unique(group))

    ᾱ ~ Normal(0, 1)
    τ ~ Exponential(1)
    α ~ filldist(Normal(ᾱ, τ), n_groups)

    σ ~ Exponential(1)

    μ = α[group]
    y ~ MvNormal(μ, σ * I)
end

simulate(α, group; σ=1.0) = [rand(Normal(α[i], σ)) for i in group]

α = [-2, 0, 2]
group = repeat(1:3, 50)
y = vec(simulate(α, group))

m = random_intercept(group, y)
chn = sample(m, NUTS(), 1000)