AR

Overview

AR is the autoregressive process of order p, whose current value is a linear combination of its p previous values plus Gaussian noise. It is the building block of latent autoregressive models, in which the states, the coefficients θ and the noise precision γ are all inferred, each with a prior of its own: a normal on θ, a gamma on γ. Chained over time, x[t] ~ AR(x[t - 1], θ, γ) in an RxInfer model, it is a state-space model whose transition is learnt.

The node is a stochastic node, a density over its variables. Its rules are variational: they compute messages from the marginals of the other variables, as variational message passing does.

Definition

For the state xₜ = (yₜ₋₁, …, yₜ₋ₚ):

\[p(y_t \mid x_t, \theta, \gamma) = \mathcal{N}\big(y_t \mid \theta^\top x_t, \, \gamma^{-1}\big).\]

For p > 1 the edges carry whole states: y = (yₜ, …, yₜ₋ₚ₊₁), related to x by the companion matrix A(θ), whose first row is θᵀ and whose sub-diagonal is ones, with noise on the first component only:

\[y = A(\theta)\, x + \varepsilon, \qquad A(\theta) = \begin{pmatrix} \theta_1 & \theta_2 & \cdots & \theta_{p-1} & \theta_p \\ 1 & 0 & \cdots & 0 & 0 \\ & \ddots & & & \vdots \\ 0 & 0 & \cdots & 1 & 0 \end{pmatrix}, \qquad \varepsilon \sim \mathcal{N}\big(0, \operatorname{diag}(\gamma^{-1}, 0, \dots, 0)\big).\]

The other components of y are copies of those of x, so the node's density is degenerate in them; the rules and the average energy account for that.

Interfaces

using AutoregressiveMessagePassingRules, MessagePassingRulesBase, BayesBase, ExponentialFamily
MessagePassingRulesBase.nodespec(AR)

The node has four interfaces; y also answers to out.

Under ARVMP(Univariate, 1, …) every edge is a scalar; under ARVMP(Multivariate, p, …) y, x and θ are vectors of length p, even for p = 1. Messages of the other form raise a MethodError. The joint marginal q(y, x) is a multivariate normal of dimension 2p, y first.

Algorithm

The node's rules run under the algorithm ARVMP(form, order, stype), which the model must name: the node declares no algorithm of its own, and under DefaultAlgorithm it has no rule. All three arguments are required, without defaults:

• form: Univariate for an AR(1) with scalar edges, or Multivariate;
• order: the order p; Univariate forces it to 1, with a warning for any other value;
• stype: ARsafe(), which forms the joint q(y, x) from its precision, regularising the noiseless components with a large finite precision, or ARunsafe(), which forms it from its covariance, exactly, through the Kalman gain. The two agree up to that regularisation, and exactly for a univariate AR(1).

ARVMP declares no dependencies: dependencies_spec returns nothing, and each rule takes the inputs of the default scheme.

Supported rules

MessagePassingRulesBase.rule_coverage(AR)

Every rule is under ARVMP, none under the default algorithm. Each target has two rules, one for each factorisation:

• the structured q(y, x) q(θ) q(γ), in which y and x exchange messages by belief propagation and the rules towards θ and γ read the joint q(y, x);
• the mean field q(y) q(x) q(θ) q(γ), in which every rule reads marginals only.

The joint marginal q(y, x) has a rule from the messages on y and x. Both factorisations have an average energy, so the node contributes to the Bethe free energy.

Example

The message towards γ of an AR(2), from the joint marginal of the states and the marginal of the coefficients: a gamma of shape 3/2 whose rate is half the expected squared residual.

@call_message_update_rule(
    node = AR, target = :γ, algorithm = ARVMP(Multivariate, 2, ARsafe()),
    clusters = ((:y, :x) => MvNormalMeanCovariance(ones(4), [1.0 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]),),
    q = (θ = MvNormalMeanCovariance([0.5, 0.25], [1.0 0.0; 0.0 1.0]),),
)

(algo, args)->begin
        order = getorder(algo)
        F = getvform(algo)
        (y_x_mean, y_x_cov) = mean_cov(args.q[:y, :x])
        (mθ, Vθ) = mean_cov(args.q[:θ])
        mA = as_companion_matrix(mθ)
        (my, Vy) = (ar_slice(F, y_x_mean, 1:order), ar_slice(F, y_x_cov, 1:order, 1:order))
        (mx, Vx) = (ar_slice(F, y_x_mean, order + 1:2order), ar_slice(F, y_x_cov, order + 1:2order, order + 1:2order))
        Vxy = ar_slice(F, y_x_cov, order + 1:2order, 1:order)
        C = add_outer(Vx, mx)
        R = add_outer(Vy, my)
        L = gaussian_cross_moment(Vxy, mx, my)
        B = (first(R) - 2 * first(mA * L)) + first(mA * C * mA') + trace_product(Vθ, C)
        GammaShapeRate(convert(eltype(B), 3 // 2), B / 2)
    end

@call_message_update_rule runs the rule and draws the marginals it read, the joint q(y, x) among them. getresult extracts the message:

julia> result = @call_message_update_rule(
           node = AR, target = :γ, algorithm = ARVMP(Multivariate, 2, ARsafe()),
           clusters = ((:y, :x) => MvNormalMeanCovariance(ones(4), [1.0 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]),),
           q = (θ = MvNormalMeanCovariance([0.5, 0.25], [1.0 0.0; 0.0 1.0]),),
       );

julia> shape(getresult(result)) ≈ 1.5 && rate(getresult(result)) ≈ 2.6875
true

Limitations

• A model must give ARVMP; without it the node has no rule.
• The rules are variational only: there is no belief propagation towards θ or γ, so both need an initial marginal, as does any edge a mean-field schedule reads first.
• The messages on y and x must be normals of the algorithm's variate form.

API

The node and its alias:

AR

julia> algorithm = ARVMP(Univariate, 1, ARsafe());

julia> result = @call_message_update_rule(
           node = AR, target = :y, algorithm = algorithm,
           m = (x = NormalMeanVariance(1.0, 1.0),),
           q = (θ = NormalMeanVariance(1.0, 1.0), γ = GammaShapeRate(1.0, 1.0)),
       );

julia> mean_var(getresult(result)) .≈ (0.5, 1.5)
(true, true)

Autoregressive

Its algorithm, and the two ways it forms the joint q(y, x):

ARVMP(form, order, stype)

julia> ARVMP(Multivariate, 3, ARsafe()) isa ARVMP
true

ARsafe()

ARunsafe()

The vector that picks one component, which the rules use for the first component of the state:

AutoregressiveMessagePassingRules.StandardBasisVector(length::Int, index::Int, scale = 1)

julia> e = AutoregressiveMessagePassingRules.StandardBasisVector(3, 2)
3-element AutoregressiveMessagePassingRules.StandardBasisVector{Int64}:
 0
 1
 0

julia> using LinearAlgebra: dot

julia> dot(e, [10.0, 20.0, 30.0])
20.0