GitHub

Getting started

This page builds a small model by hand and runs inference on it with the engine alone. It does what RxInfer's @model and infer do for you. The model is a latent x with a normal prior, observed through normal noise of known variance:

\[x \sim \mathcal{N}(0, 10), \qquad y \mid x \sim \mathcal{N}(x, 1).\]

Its posterior is normal, with precision $1/10 + 1$ and mean $y / (1/10 + 1)$. Its evidence is $p(y) = \mathcal{N}(y \mid 0, 11)$. Belief propagation computes both exactly, and the Bethe free energy of the graph equals $-\log p(y)$.

Packages and the node

The engine knows no factor node until one is declared. Nodes and their rules come from rule packages, which declare them with MessagePassingRulesBase. This site declares one small node of its own, Gaussian, a normal density with a known variance, in a file that every page includes (The example node shows it). Rocket.jl provides subscribe!.

using ReactiveMP, MessagePassingRulesBase, BayesBase, ExponentialFamily, Distributions, Rocket
import ReactiveMP: activate!, FactorNodeActivationOptions, get_stream_of_marginals

include(joinpath(pkgdir(ReactiveMP), "docs", "nodes.jl"))

The declaration names the node's interfaces: out, the output, μ, the mean, and v, the variance.

MessagePassingRulesBase.nodespec(Gaussian)
NodeSpec: Gaussian
μvoutGaussianstochastic
interfacealiases
out
μ
v
default algorithmMessagePassingRulesBase.DefaultAlgorithm()
static inputsnone
defineddocs/nodes.jl:11

Both factors of the model are this node: the prior is $\mathcal{N}(x \mid 0, 10)$ and the likelihood $\mathcal{N}(y \mid x, 1)$.

1. Variables

A randomvar is inferred. A datavar receives observations. A constvar holds a fixed value, such as a known mean or variance (see Variables).

x = randomvar(label = :x)
y = datavar(label = :y)
prior_mean, prior_var, noise_var = constvar(0.0), constvar(10.0), constvar(1.0)

2. Factor nodes

A factornode connects a node type to its variables, one (interface, variable) pair per edge (see Factor nodes).

prior = factornode(Gaussian, [(:out, x), (:μ, prior_mean), (:v, prior_var)])
likelihood = factornode(Gaussian, [(:out, y), (:μ, x), (:v, noise_var)])
FactorNode Gaussian (stochastic)
  out ── y    data
  μ   ── x    random
  v   ── 1.0  constant
  clusters: (out, μ, v)

The likelihood lists its edges, the variable on each and its kind, and its clusters. Without a factorisation, all of a node's interfaces form one cluster, (out, μ, v), and the engine runs the node's belief propagation rules.

3. Activation

Activation wires the streams, variables first, then nodes (see Inference lifecycle). Each takes its options: RandomVariableActivationOptions, DataVariableActivationOptions and ReactiveMP.FactorNodeActivationOptions. Here the nodes track log scales (see Activation options). A constant needs no activation.

activate!(x, RandomVariableActivationOptions())
activate!(y, DataVariableActivationOptions())
for node in (prior, likelihood)
    activate!(node, FactorNodeActivationOptions(; logscales = true))
end

4. Subscribing

The engine computes nothing until something listens. Subscribe to the posterior of x, and to the free energy. bethe_free_energy takes every node and every variable of the graph, the constants included.

posteriors = Marginal[]
energies = Float64[]
variables = (x, y, prior_mean, prior_var, noise_var)
subscriptions = [
    subscribe!(get_stream_of_marginals(x), (q) -> push!(posteriors, q)),
    subscribe!(bethe_free_energy(Float64, (prior, likelihood), variables), (f) -> push!(energies, f)),
]

5. Observing

new_observation! gives y a value, which propagates through the graph. The likelihood computes its message towards x. The variable x multiplies it with the prior's message, and the posterior and the free energy update.

new_observation!(y, 2.0)
q = last(posteriors)
Marginal(ExponentialFamily.NormalWeightedMeanPrecision{Float64}(xi=2.0, w=1.1)) with logscale = -2.2997043514220397

The posterior is in its natural form: weighted mean $y = 2$ and precision $1.1$. It matches the closed form, mean $2 / 1.1$ and variance $1 / 1.1$:

mean(q) ≈ 2.0 / 1.1, var(q) ≈ 1 / 1.1
(true, true)

The posterior's log scale (see Log scales) is the log evidence, and the free energy is minus the log evidence:

log_evidence = logpdf(Normal(0, sqrt(11)), 2.0)
getlogscale(q) ≈ log_evidence, last(energies) ≈ -log_evidence
(true, true)

Each new observation updates them again. This is how the engine serves streaming data:

new_observation!(y, 1.0)
mean(last(posteriors)), length(energies)
(0.9090909090909091, 2)

When the run is over, unsubscribe:

foreach(unsubscribe!, subscriptions)

Next