ConjugateAR

Overview

ConjugateAR is the likelihood of AR with the coefficients θ and the noise precision γ joint on one edge w = (θ, γ). With a normal-gamma prior on w the parameter update is conjugate: the posterior q(w) is the normal-gamma posterior of a Bayesian linear regression with unknown noise precision, with the states' expected sufficient statistics in place of the data. Use it where θ and γ should stay coupled, rather than factorised as q(θ) q(γ).

Here q(w) is the marginal of w, the approximate posterior that message passing computes. As for AR, the node's rules are those of variational message passing.

Definition

For the state xₜ = (yₜ₋₁, …, yₜ₋ₚ) and w = (θ, γ):

\[p(y \mid x, w) = \mathcal{N}\big(y_1 \mid \theta^\top x, \, \gamma^{-1}\big),\]

with y = A(θ) x + ε in companion form, noise on the first component only, as for AR. Averaged over q(y, x), the node is a normal-gamma factor in w with statistics

\[C = \langle x x^\top \rangle, \qquad b = \langle x\, y_1 \rangle, \qquad a = \langle y_1^2 \rangle,\]

that is Λ = C, μ = C⁻¹ b, shape (3 - p)/2 and rate (a - bᵀ C⁻¹ b)/2. These are the parameters of the message towards w.

Interfaces

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

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

The joint marginal q(y, x) is a multivariate normal of dimension 2p, y first.

Algorithm

The node shares AR's algorithm, ARVMP(form, order, stype), which the model must name, with every argument required. form is Univariate for an AR(1) on scalars, whose rules read q(w)'s one coefficient as a scalar θ, or Multivariate for an AR of any order on state vectors. stype is ARsafe() or ARunsafe(), as for AR. The node declares no dependencies, so each rule takes the inputs of the default scheme.

Supported rules

MessagePassingRulesBase.rule_coverage(ConjugateAR)

Every rule is under ARVMP. The rules towards y and x exist for the structured factorisation q(y, x) q(w), as belief propagation between them, and for the mean field. They and the joint marginal q(y, x) are AR's, computed from the marginals of θ and γ that q(w) implies. The rule towards w and the average energy need the structured factorisation q(y, x) q(w).

Example

The message towards w of an AR(1) from the joint marginal of the states:

@call_message_update_rule(
    node = ConjugateAR, target = :w, algorithm = ARVMP(Multivariate, 1, ARsafe()),
    clusters = ((:y, :x) => MvNormalMeanCovariance([1.0, 0.5], [1.0 0.0; 0.0 1.0]),),
)

(algo, args)->begin
        order = getorder(algo)
        (myx, Vyx) = mean_cov(args.q[:y, :x])
        x_range = order + 1:2order
        mx = myx[x_range]
        my1 = first(myx)
        Vx = Vyx[x_range, x_range]
        Vy1 = Vyx[1, 1]
        cov_x_y1 = Vyx[x_range, 1]
        C = Vx + mx * transpose(mx)
        b = cov_x_y1 + mx * my1
        a = Vy1 + my1 ^ 2
        μ = cholinv(C) * b
        α = convert(eltype(μ), (3 - order) / 2)
        β = (a - dot(b, μ)) / 2
        MvNormalGamma(μ, C, α, β)
    end

@call_message_update_rule runs the rule and draws the joint marginal it read. With C = 1.25, b = 0.5 and a = 2, the getresult of the call is the normal-gamma of the definition:

julia> result = @call_message_update_rule(
           node = ConjugateAR, target = :w, algorithm = ARVMP(Multivariate, 1, ARsafe()),
           clusters = ((:y, :x) => MvNormalMeanCovariance([1.0, 0.5], [1.0 0.0; 0.0 1.0]),),
       );

julia> μ, Λ, α, β = params(getresult(result));

julia> μ ≈ [0.4] && Λ ≈ fill(1.25, 1, 1) && α ≈ 1.0 && β ≈ 0.9
true

Limitations

• A model must give ARVMP; without it the node has no rule.
• The rule towards w and the average energy need the structured factorisation q(y, x) q(w): under the mean field q(y) q(x) q(w) there is no message towards w and no average energy.
• The message towards w is improper for p ≥ 3, its shape (3 - p)/2 being non-positive; its product with a proper normal-gamma prior is proper. q(w) therefore needs a prior, and an initial marginal where the schedule reads it first.

API

ConjugateAR

julia> result = @call_message_update_rule(
           node = ConjugateAR, target = :w, algorithm = ARVMP(Multivariate, 1, ARsafe()),
           clusters = ((:y, :x) => MvNormalMeanCovariance([1.0, 0.5], [1.0 0.0; 0.0 1.0]),),
       );

julia> getresult(result) isa MvNormalGamma
true