PolyaMessagePassingRules

Licence: GPL-3

This package is licensed under GPL-3, because its dependency PolyaGammaHybridSamplers is. The engine and the other rule packages are MIT and do not depend on it; a model that loads this package is subject to the GPL.

Two factor nodes for count data through the logistic, BinomialPolya and MultinomialPolya. Both are augmented with Pólya-Gamma variables, so that the messages towards their weights are normal and a normal prior on the weights stays conjugate. Use them for binomial or logistic regression, and for multinomial regression.

PolyaMessagePassingRules — Module

PolyaMessagePassingRules

The Pólya-Gamma augmented nodes, for count data through the logistic:

BinomialPolya, y ~ Binomial(n, σ(xᵀβ)), a binomial (or logistic) regression with covariates x and weights β, under BinomialPolyaApproximation;
MultinomialPolya, x ~ Multinomial(N, p(ψ)) through logistic stick-breaking, under MultinomialPolyaApproximation, with the helpers logistic_stick_breaking and compose_Nks.

Pólya-Gamma augmentation makes the message towards the weights normal, so a normal prior on them stays conjugate. Each rule towards the weights reads the message on its own edge, which a model must initialise. Both nodes have average energies.

This package is licensed under GPL-3, through its dependency PolyaGammaHybridSamplers, which is. The engine and the other rule packages do not depend on it.

Examples

julia> result = @call_message_update_rule(
           node = BinomialPolya, target = :β,
           q = (y = PointMass(3), x = PointMass(0.1), n = PointMass(5)),
           m = (β = NormalWeightedMeanPrecision(0.0, 1.0),),
       );

julia> weightedmean(getresult(result)) ≈ 0.05 && precision(getresult(result)) ≈ 0.0125
true

The weighted mean is (y - n / 2) x and the precision ω x², ω = n / 4 being the Pólya-Gamma mean at xᵀβ = 0.

source

Where these rules run

ReactiveMP runs these rules on a factor graph, which RxInfer builds from a GraphPPL model. The examples here call the rules directly, as a test does.

The augmentation

For a count y out of b trials with log-odds ψ, the Pólya-Gamma identity of Polson, Scott and Windle (2013),

\[\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 ω ~ PG(b, ψ). The rules replace ω by its mean at the current estimate of ψ, read from the message on the weights' own edge. A model must therefore give that message an initial message. The average energies are those of the original likelihoods, not of the augmented ones.

PolyaMessagePassingRules

The augmentation

Pages