The Probability Interface

The easiest way to manipulate and query DynamicPPL models is via the DynamicPPL probability interface.

Let's use a simple model of normally-distributed data as an example.

using DynamicPPL
using Distributions
using FillArrays

using LinearAlgebra
using Random

@model function gdemo(n)
    μ ~ Normal(0, 1)
    x ~ MvNormal(Fill(μ, n), I)
    return nothing
end

We generate some data using μ = 0:

Random.seed!(1776)
dataset = randn(100)

Conditioning and Deconditioning

Bayesian models can be transformed with two main operations, conditioning and deconditioning (also known as marginalization). Conditioning takes a variable and fixes its value as known. We do this by passing a model and a collection of conditioned variables to | or its alias condition:

model = gdemo(length(dataset)) | (x=dataset, μ=0)

This operation can be reversed by applying decondition:

decondition(model)

We can also decondition only some of the variables:

decondition(model, :μ)

gdemo(x::AbstractVector{<:Real}) = gdemo(length(x)) | (; x)

Probabilities and Densities

We often want to calculate the (unnormalized) probability density for an event. This probability might be a prior, a likelihood, or a posterior (joint) density. DynamicPPL provides convenient functions for this. For example, we can calculate the joint probability of a set of samples (here drawn from the prior) with logjoint:

model = gdemo(length(dataset)) | (x=dataset,)

Random.seed!(124)
sample = rand(model)
logjoint(model, sample)

-181.7247437162069

For models with many variables rand(model) can be prohibitively slow since it returns a NamedTuple of samples from the prior distribution of the unconditioned variables. We recommend working with samples of type DataStructures.OrderedDict in this case:

using DataStructures

Random.seed!(124)
sample_dict = rand(OrderedDict, model)
logjoint(model, sample_dict)

-181.7247437162069

The prior probability and the likelihood of a set of samples can be calculated with the functions loglikelihood and logjoint, respectively:

logjoint(model, sample) ≈ loglikelihood(model, sample) + logprior(model, sample)

true

logjoint(model, sample_dict) ≈
loglikelihood(model, sample_dict) + logprior(model, sample_dict)

true

Example: Cross-validation

To give an example of the probability interface in use, we can use it to estimate the performance of our model using cross-validation. In cross-validation, we split the dataset into several equal parts. Then, we choose one of these sets to serve as the validation set. Here, we measure fit using the cross entropy (Bayes loss).^[1] (For the sake of simplicity, in the following code, we enforce that nfolds must divide the number of data points. For a more competent implementation, see MLUtils.jl.)

# Calculate the train/validation splits across `nfolds` partitions, assume `length(dataset)` divides `nfolds`
function kfolds(dataset::Array{<:Real}, nfolds::Int)
    fold_size, remaining = divrem(length(dataset), nfolds)
    if remaining != 0
        error("The number of folds must divide the number of data points.")
    end
    first_idx = firstindex(dataset)
    last_idx = lastindex(dataset)
    splits = map(0:(nfolds - 1)) do i
        start_idx = first_idx + i * fold_size
        end_idx = start_idx + fold_size
        train_set_indices = [first_idx:(start_idx - 1); end_idx:last_idx]
        return (view(dataset, train_set_indices), view(dataset, start_idx:(end_idx - 1)))
    end
    return splits
end

function cross_val(
    dataset::Vector{<:Real};
    nfolds::Int=5,
    nsamples::Int=1_000,
    rng::Random.AbstractRNG=Random.default_rng(),
)
    # Initialize `loss` in a way such that the loop below does not change its type
    model = gdemo(1) | (x=[first(dataset)],)
    loss = zero(logjoint(model, rand(rng, model)))

    for (train, validation) in kfolds(dataset, nfolds)
        # First, we train the model on the training set, i.e., we obtain samples from the posterior.
        # For normally-distributed data, the posterior can be computed in closed form.
        # For general models, however, typically samples will be generated using MCMC with Turing.
        posterior = Normal(mean(train), 1)
        samples = rand(rng, posterior, nsamples)

        # Evaluation on the validation set.
        validation_model = gdemo(length(validation)) | (x=validation,)
        loss += sum(samples) do sample
            logjoint(validation_model, (μ=sample,))
        end
    end

    return loss
end

cross_val(dataset)

-212760.30282411768