BinomialPolya

BinomialPolyaApproximation(; samples = nothing, points = 32)

The algorithm of BinomialPolya, and its default: a model names it only to sample or to change points.

Keywords

• samples: nothing or a number of draws. Default nothing. With nothing, the rules use means: the Pólya-Gamma mean at xᵀ times the mean of β, and σ(xᵀβ) at the mean of β. With a number k, they average over k draws of β, the rule towards β over one Pólya-Gamma draw for each, from the rule context's generator, ctx.rng. An engine supplies it; a call by hand passes ctx = MessagePassingRulesBase.RuleContext(rng = …).
• points: the number of Gauss–Hermite points with which the average energy computes ⟨softplus(ψ)⟩. Default 32. The messages do not depend on it.

The average energy does not sample under either: xᵀβ is normal under a normal q(β), so its expectation is computed by Gauss–Hermite cubature.

Examples

julia> BinomialPolyaApproximation().samples === nothing
true

julia> BinomialPolyaApproximation(samples = 100).samples
100

julia> BinomialPolyaApproximation(points = 64).points
64

BinomialPolya

The stochastic node of a binomial regression through the logistic,

\[p(y \mid x, n, \beta) = \binom{n}{y} \sigma(x^\top \beta)^{y} \big(1 - \sigma(x^\top \beta)\big)^{n - y}, \qquad \sigma(\psi) = \frac{1}{1 + e^{-\psi}},\]

with its weights β Pólya-Gamma augmented so that their message is normal. With n = 1 it is logistic regression.

Interfaces

• y: the number of successes, observed, a PointMass;
• x: the covariates, observed, a PointMass holding a vector (or a number, for a scalar β);
• n: the number of trials, observed, a PointMass;
• β: the weights, with a normal message and marginal, MvNormal for a vector x.

Augmentation

For a count y out of b trials with log-odds ψ, the Pólya-Gamma identity

\[\frac{(e^{\psi})^{y}}{(1 + e^{\psi})^{b}} = 2^{-b} e^{\kappa \psi} \int_0^\infty e^{-\omega \psi^2 / 2} \, \mathrm{PG}(\omega \mid b, 0) \, d\omega, \qquad \kappa = y - b / 2,\]

makes the likelihood Gaussian in ψ given the auxiliary ω ~ PG(b, ψ): proportional to exp(κψ - ωψ²/2). The rules replace ω by its expectation, E[ω] = b / (2c) · tanh(c / 2) for PG(b, c), at the current estimate c of ψ, so the message towards the weights is a normal in weighted-mean and precision form.

With ψ = xᵀβ and b = n, the message towards β has the weighted mean (y - n/2) x and the precision ω x xᵀ, where ω is the Pólya-Gamma mean at xᵀ times the mean of the message on β. The message towards y is Binomial(n, σ(xᵀ mᵦ)), mᵦ the mean of q(β). Both rules average over draws of β instead when BinomialPolyaApproximation is given a number of samples.

The average energy is that of the binomial likelihood itself, not of the augmented one: -log C(n, y) - y ⟨ψ⟩ + n ⟨softplus(ψ)⟩, with ⟨softplus(ψ)⟩ under the normal ψ = xᵀβ by Gauss–Hermite cubature with the algorithm's points.

Limitations

• The rule towards β reads the message on β itself, besides the marginals its factorisation gives it, so a model must initialise that message.
• y, x and n must be observed: their rules and the average energy take PointMass marginals only, and there is no rule towards x or n.
• The Gauss–Hermite cubature of ⟨softplus(ψ)⟩ is exact to rounding while ψ is narrow, a variance up to about 1, and loses accuracy as ψ broadens, since softplus bends sharply on the cubature's scale: about 1e-4 relative at a variance of 25 and 1e-2 at 400 with 32 points. More points help slowly. The free energy of a model with broad priors is therefore approximate in its first iterations; the messages do not use the cubature.

Examples

julia> result = @call_message_update_rule(
           node = BinomialPolya, target = :y,
           q = (x = PointMass([1.0, 2.0]), n = PointMass(10), β = MvNormalMeanCovariance([0.0, 0.0], [1.0 0.0; 0.0 1.0])),
       );

julia> mean(getresult(result)) ≈ 5.0
true

At β = 0 the success probability is σ(0) = 1/2, so the message is Binomial(10, 1/2).

See also BinomialPolyaApproximation, MultinomialPolya.