Guides / Choose a prior

Choose a prior¶

The prior on the population matters most when the data say least: many bins, few objects, or measurement errors as wide as the bins. Start with the default. Change it when the checks at the end of this page say the prior, not the data, sets the answer. What each prior means is in the theory.

Your situation	Prior to try
no strong reason to choose	`Dirichlet(1.0)`, the default
about as many bins as objects	`Dirichlet.reference()`
a jagged fit, with errors as wide as the bins	`LogisticGP()` or `LogisticRandomWalk(order=2)`
a reference sample of the same population	`Dirichlet.from_counts(counts, strength)`
unequal bins, with a Dirichlet	`Dirichlet.from_base_measure(...)`
expected numbers of objects, with prior information on them	`GammaRates(shape, rate)`

binference’s priors fall into three families. The Dirichlet prior family is conjugate, so every engine takes it. The logistic priors add smoothness across neighbouring bins; Gibbs and VB do not accept them. The Gamma rates prior models expected counts rather than fractions.

import numpy as np
import binference
from binference import priors

mock = binference.simulate.mock_1d(1000, binference.Bins.linear(0, 3, 30), noise=(0.1, 0.3), rng=4)
lk = mock.likelihoods

Pass a prior as the second argument of infer(), or as prior= to PopulationModel and PopulationEstimator.

The Dirichlet family¶

flat = priors.Dirichlet(1.0)                         # the default
sparse = priors.Dirichlet.reference()                # concentration 1/K per bin
by_volume = priors.Dirichlet.from_base_measure(5.0)  # concentration in proportion to bin width
ref_counts = np.bincount(mock.truth.labels[:50], minlength=30)
informed = priors.Dirichlet.from_counts(ref_counts, strength=0.5)

for p in (flat, sparse, informed):
    print(p, binference.infer(lk, p, method="map").rho[:3].round(4))

Dirichlet(1.0): use it unless you have a reason not to.
Dirichlet.reference(): use it when the number of bins \(K\) is comparable to the number of objects \(N\). There the default biases the fit toward equal fractions.
Dirichlet.from_base_measure(a): use it on unequal bins, so the prior means the same density however you bin.
Dirichlet.from_counts(counts, strength): use it when you have the counts \(n^{\rm ref}_k\) of a reference sample. strength is \(\kappa\): 1 trusts the sample fully, and smaller values trust it less.

GeneralizedDirichlet, DirichletTree and PolyaTree stay conjugate, with more freedom in how bins co-vary:

gd = priors.GeneralizedDirichlet(a=2.0, b=6.0)        # stick-breaking Betas
tree = priors.DirichletTree([[0, 1, 2], [3, [4, 5]]], alpha=2.0)
polya = priors.PolyaTree(c=2.0)                       # dyadic, multi-scale
print(binference.infer(lk, polya, method="map").rho[:3].round(4))

DirichletTree: use it when groups of bins belong together, such as subclasses of a class. Pass the groups as nested lists of bin indices, or as dicts of named groups; on categorical bins a bin can be named by its label. A dict alpha sets the concentration of the branch into each group, or bin as the tree writes it, that it names; the other branches keep 1.
PolyaTree: use it on 1-D ordered bins for a conjugate prior with structure at several scales.

describe(bins) draws the tree of any prior in this family: each branch with its concentration, the share of its parent it takes on average, and the mean fraction of the population below it. Pass counts to see the posterior given them.

classes = binference.Bins.categories(["Sa", "Sb", "E0", "E3", "S0"])
early = {"elliptical": ["E0", "E3"], "lenticular": ["S0"]}
named = priors.DirichletTree({"spiral": ["Sa", "Sb"], "early": early}, alpha={"spiral": 3.0})
print(named.describe(classes))

DirichletTree: prior, on 5 bins
root: Sa, Sb, E0, E3, S0
├── α=3  share 0.75  mass 0.750  spiral: Sa, Sb
│   ├── α=1  share 0.50  mass 0.375  Sa
│   └── α=1  share 0.50  mass 0.375  Sb
└── α=1  share 0.25  mass 0.250  early: E0, E3, S0
    ├── α=1  share 0.50  mass 0.125  elliptical: E0, E3
    │   ├── α=1  share 0.50  mass 0.062  E0
    │   └── α=1  share 0.50  mass 0.062  E3
    └── α=1  share 0.50  mass 0.125  lenticular: S0

Smoothness: the logistic family¶

If the Dirichlet fit is jagged and the errors are as wide as the bins, use a smoothness prior. It infers how smooth the population is:

rw = priors.LogisticRandomWalk(order=2)
gp = priors.LogisticGP(kernel="matern32")
print(binference.infer(lk, gp, method="map").rho[:3].round(4))

prior	use it when	hyperparameters (inferred by default)
`LogisticRandomWalk(order=1)`	neighbouring bins should be close	`scale`
`LogisticRandomWalk(order=2)`	the population follows smooth trends	`scale`
`LogisticGP(kernel=...)`	you want a stationary Gaussian process: `"matern12"`, `"matern32"`, `"matern52"`, `"rbf"`	`amplitude`, `length_scale` (one per axis)
`LogisticNormal(mean, cov)`	you have a covariance, or want independent bins	`scale`

Grids. All four work on 1-D bins and on grids. Only LogisticNormal accepts masked grids and categories.
Rough populations. If a LogisticGP fit looks too smooth, raise n_basis (defaults).

Hyperpriors. Give a hyperparameter a hyperprior, or a number to fix it:

gp_fixed = priors.LogisticGP(amplitude=1.5, length_scale=priors.LogNormal(np.log(0.5), 0.5))
rw_tight = priors.LogisticRandomWalk(order=2, scale=priors.HalfNormal(0.3))

A reference population. With reference its fraction in each bin, mean=np.log(reference) centres the prior on it.

Rates instead of fractions¶

GammaRates(shape, rate) puts independent Gamma priors on the expected number of objects per bin, \(\lambda_k\). The data then form a Poisson process. Every fit of fractions already gives the total and the rates, under a scale-free prior on the total (Selection and rates). Use GammaRates when the rates need their own prior, for example bin-specific shapes or an informative total. Combine it with an exposure or a selection function. Occurrence rates per star, or counts per unit time or area, are estimated this way.

est = binference.infer(lk, priors.GammaRates(shape=1.0, rate=1e-3), method="map")
print(est.rates.sum().round(1), est.rho.sum())

Sampled results hold both rates, the draws of \(\lambda\), and samples, the fractions they imply. NUTS and Gibbs both sample rates.

Check a prior¶

Draw from the prior to see which populations it allows before any data:

model = binference.PopulationModel(lk, gp)
draws = model.sample_prior(5, rng=0)       # five populations from the prior
print(draws.shape, np.allclose(draws.sum(axis=1), 1.0))

After a fit, two numbers say how much the prior sets the answer. A bin with little data support takes its fraction from the prior (Choose bins). result.prior_information_ratio() compares the prior’s information with the data’s. It is small when the data dominate, and near one or above when the prior does.

fit = binference.infer(lk, rw, method="map")
print(round(fit.prior_information_ratio(), 3))

To choose between priors, fit each and compare the log predictive density of held-out objects (Read a result). binference.simulate.prior_predictive() simulates whole catalogues from the prior, which is what simulation-based calibration needs (Check a result).