MultinomialPolya

MultinomialPolyaApproximation(; points = 21)

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

Keywords

• points: the number of Gauss–Hermite points with which the average energy computes each ⟨softplus(ψ_k)⟩. Default 21. The messages do not depend on it.

Examples

julia> MultinomialPolyaApproximation().points
21

MultinomialPolya

The stochastic node of a multinomial through logistic stick-breaking,

\[p(x \mid N, \psi) = \frac{N!}{\prod_{k=1}^{K} x_k!} \prod_{k=1}^{K} p_k^{x_k}, \qquad p_k = \sigma(\psi_k) \prod_{j<k} \big(1 - \sigma(\psi_j)\big), \quad p_K = \prod_{j<K} \big(1 - \sigma(\psi_j)\big),\]

for K categories and K - 1 log-odds ψ, Pólya-Gamma augmented so that the message towards ψ is normal. The multinomial factors into K - 1 binomials, the k-th of x_k out of the N_k = N - Σ_{j<k} x_j trials left, with log-odds ψ_k. With ψ a linear function of covariates it is multinomial regression.

Interfaces

• x: the counts of the K categories, a vector: observed (PointMass) or a Multinomial;
• N: the number of trials, a PointMass, or a Poisson, Binomial or Categorical, whose mode the rules use;
• ψ: the K - 1 log-odds, MvNormal, or a univariate Normal for K = 2.

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.

Applied to each binomial of the stick with b = N_k, the message towards ψ has the weighted mean x_k - N_k / 2 and the diagonal precision ω_k, the Pólya-Gamma mean at the k-th entry of the mean of the message on ψ. The message towards x is the multinomial at logistic_stick_breaking of the mean of q(ψ).

The average energy is that of the multinomial itself, -⟨log C(x)⟩ - Σ_{k<K} ⟨x_k⟩⟨ψ_k⟩ + Σ_{k<K} N_k ⟨softplus(ψ_k)⟩, with each ⟨softplus(ψ_k)⟩ by Gauss–Hermite cubature over the mean and variance of ψ_k, with the points of MultinomialPolyaApproximation. For a Multinomial q(x), the coefficient's expectation is log N! - Σ_k ⟨log x_k!⟩ over its binomial marginals.

Limitations

• The rule towards ψ reads the message on ψ itself, besides the marginals its factorisation gives it, so a model must initialise that message.
• There is no rule towards N, and the average energy takes a PointMass q(N) only.
• The rules use the mean of ψ only: its variance enters the average energy, not the messages.
• 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 = MultinomialPolya, target = :ψ,
           q = (x = PointMass([2, 3, 5]), N = PointMass(10)),
           m = (ψ = MvNormalWeightedMeanPrecision(zeros(2), [1.0 0.0; 0.0 1.0]),),
       );

julia> weightedmean(getresult(result)) ≈ [-3.0, -1.0]
true

julia> precision(getresult(result)) ≈ [2.5 0.0; 0.0 2.0]
true

The stick's two breaks leave N_k = 10 and 8 trials, so the weighted means are 2 - 5 and 3 - 4 and the precisions, at ψ = 0, N_k / 4.

See also MultinomialPolyaApproximation, logistic_stick_breaking, compose_Nks, BinomialPolya.

logistic_stick_breaking(m) -> Vector{Float64}

The K probabilities of logistic stick-breaking with the K - 1 log-odds m: p_k = σ(m_k) Π_{j<k} (1 - σ(m_j)) for k < K, and the rest of the stick, Π_{j<K} (1 - σ(m_j)), for the last. The probabilities of MultinomialPolya at ψ = m.

Examples

julia> using PolyaMessagePassingRules

julia> logistic_stick_breaking([0.0, 0.0])
3-element Vector{Float64}:
 0.5
 0.25
 0.25

compose_Nks(x, N) -> Vector

The number of trials left at each of the first K - 1 breaks of the stick, for the K counts x out of N: N_k = N - Σ_{j<k} x_j, of the element type of x. The trials of the binomials MultinomialPolya factors into.

Examples

julia> using PolyaMessagePassingRules

julia> compose_Nks([2, 3, 5], 10)
2-element Vector{Int64}:
 10
  8