Pólya Trees

Python · R  ·  Download Pólya-tree module

A Prior on the Distribution Itself

The arc closes on the most ambitious prior in it: one over whole distributions. A Pólya tree splits the real line in half, then splits each half, and so on, placing a Beta prior on the probability of going left at every node. Centre those Betas on a parametric family and the prior sits on that family; make them tight and the posterior can barely move; make them loose and the data can drag the density anywhere. The single parameter cc that controls the tightness is, quite literally, how much you trust the parametric model.

Yε0∼Beta(c m2, c m2),G(Bε0)=G(Bε) Yε0at level mY_{\varepsilon 0}\sim\text{Beta}(c\,m^2,\,c\,m^2),\qquad G(B_{\varepsilon 0})=G(B_\varepsilon)\,Y_{\varepsilon 0}\quad\text{at level }m

As a density estimator

That construction pays twice. As a density estimator it starts from a normal centre and lets the data override it: on a bimodal mixture the tree recovers both modes with an L1L_1 error of 0.191 against a kernel estimate's 0.309, and unlike the KDE it returns a credible band. The centring normal it started from is drawn on the same axes — and is exactly the null hypothesis the next section tests.

As a test of the family at its centre

Because the Pólya tree is conjugate, the marginal likelihood is available in closed form, which turns the same object into a goodness-of-fit test: a Bayes factor against the parametric family at its centre. Validated on three simulated samples, it agrees with the classical tests where they are confident — normal data passes at log⁡BF=−15.5\log\text{BF}=-15.5, gamma-skewed and t3t_3 data are rejected decisively — but it does something the tests cannot. A rejected KS test tells you the normal is wrong. The Pólya tree hands back the density that fits instead.

simulated samplePólya-tree log⁡BF\log\text{BF}verdictKS ppShapiro pp
truly normal−15.5normal OK0.930.56
skewed (gamma)88.8reject5.0e−082.4e−23
heavy-tailed t3t_328.9reject6.6e−055.5e−21
S&P 500 daily returns257reject5.1e−333.7e−44

The Closing Argument on Returns

Turned on 3,772 daily S&P 500 returns, it delivers the arc's closing argument. The log Bayes factor against normality is 257. That is the log of the ratio of the two models’ marginal likelihoods, on a scale where 5 already counts as decisive, so this is not a marginal rejection but an annihilation — with skew −0.72 and excess kurtosis 11.5. The fitted density has a taller, sharper peak and tails sitting orders of magnitude above the Gaussian on a log scale. That is precisely the leptokurtosis the Risk and Asset Allocation arc spends eight projects modelling with Student-tt GARCH, extreme-value theory and copulas; here a prior over distributions confirms and quantifies it without assuming any of them.

Which method had to assume continuity

One assumption got checked rather than assumed, because R warned about it and the notebooks had passed over the warning. KS and Shapiro both assume a continuous distribution — no repeated values — and this series has 11 repeated observations out of 3,772, every one an exact zero: days the index closed unchanged. At p=5×10−33p=5\times10^{-33} the verdict is nowhere near sensitive to eleven ties, so nothing here changes. But it is worth noting which method needed the assumption. The Pólya tree partitions the line and counts what lands in each cell, so repeated values violate nothing it relies on — a small structural advantage, and one that matters on price data, where exact zeros and rounded ticks are ordinary rather than pathological.

Notebooks

Python Densities drawn from the prior at several trust levels, density estimation against a kernel estimate on a bimodal mixture, the Bayes-factor normality test validated on three known samples, and the S&P 500 verdict (Python) Open → R Independent from-scratch base-R implementation — the same prior draws, density estimate, Bayes factor and comparison against ks.test and shapiro.test (R) Open →

Downloads

bnp_polyatree.py Pólya-tree fitting by dyadic partition of the centring CDF, posterior density with credible bands, and the closed-form marginal-likelihood Bayes factor against the centring family (NumPy / SciPy) spy_returns.csv S&P 500 daily returns — 3,772 trading days, 2010–2024

Pólya-Tree Module — Source Code

"""
polyatree.py -- POLYA TREES: a nonparametric prior on a distribution (from scratch).

Backs the notebooks in  "Polya Trees".

A Dirichlet process is a prior over DISCRETE distributions. A Polya tree (Lavine 1992;
Ferguson 1974) is a prior over CONTINUOUS ones -- a random density -- built by recursively
splitting the sample space and randomising how mass flows down the splits.

Recipe. Map the data to [0,1] through the CDF of a centring distribution G0, u_i = G0(x_i).
Split [0,1] dyadically: level m has 2^m intervals, each parent splitting into two children.
At every split, the fraction of the parent's mass going LEFT is a Beta random variable

    Y_eps ~ Beta(alpha_m, alpha_m),   alpha_m = c * m^2   (Lavine's canonical choice),

independent across splits. alpha_m = c m^2 makes the random measure absolutely continuous (a
genuine density), and the concentration c controls how tightly it hugs G0: c -> infinity gives
back G0 exactly, small c allows wild departures. Because the splits are at G0's dyadic quantiles
(uniform on the u-scale), the prior is CENTRED on G0: E[density] = G0.

Conjugacy makes it trivial. With counts n_{eps0}, n_{eps1} of observations falling in the two
children of a node, the posterior is another Polya tree:

    Y_eps | data ~ Beta(alpha_m + n_{eps0}, alpha_m + n_{eps1}).

So the posterior density is available in closed form -- no MCMC, just Beta draws for the bands.

GOODNESS OF FIT. Centring on a fitted parametric G0 turns the Polya tree into a test of that
family: if G0 is right the data are uniform on the u-scale (every Y ~ 1/2) and the tree adds
nothing; if not, the branch counts pull the Y's away from 1/2. The Bayes factor against G0 is a
closed-form sum over nodes -- large values are evidence the parametric family is wrong.
"""

import numpy as np
from scipy.special import betaln


# --------------------------------------------------------------------------- #
#  fit (store dyadic counts of the CDF-transformed data)                        #
# --------------------------------------------------------------------------- #

def polyatree_fit(x, cdf, pdf, c=1.0, M=8):
    """Centre a Polya tree of depth M on the distribution with the given cdf/pdf, concentration c.
    Stores the counts of the transformed data in every dyadic interval up to level M."""
    x = np.asarray(x, float)
    u = np.clip(cdf(x), 1e-12, 1 - 1e-12)
    counts = [np.bincount(np.minimum((u * 2 ** m).astype(int), 2 ** m - 1),
                          minlength=2 ** m).astype(float) for m in range(M + 1)]
    return dict(counts=counts, c=c, M=M, cdf=cdf, pdf=pdf)


def _leaf_masses(pt, draw=False, rng=None):
    """mass in each of the 2^M leaf intervals (posterior mean, or a posterior draw)."""
    M, c = pt["M"], pt["c"]; mass = np.ones(1)
    for m in range(M):
        a = c * (m + 1) ** 2; nc = pt["counts"][m + 1]
        n0 = nc[0::2]; n1 = nc[1::2]
        y = rng.beta(a + n0, a + n1) if draw else (a + n0) / (2 * a + n0 + n1)
        nm = np.empty(2 ** (m + 1)); nm[0::2] = mass * y; nm[1::2] = mass * (1 - y)
        mass = nm
    return mass


# --------------------------------------------------------------------------- #
#  posterior density on a grid                                                  #
# --------------------------------------------------------------------------- #

def polyatree_density(pt, xgrid, ndraws=0, rng=None):
    """posterior density on xgrid: mean, and (if ndraws>0) an (ndraws, len(grid)) sample."""
    M = pt["M"]; xgrid = np.asarray(xgrid, float)
    u = np.clip(pt["cdf"](xgrid), 1e-12, 1 - 1e-12)
    leaf = np.minimum((u * 2 ** M).astype(int), 2 ** M - 1)
    g0 = pt["pdf"](xgrid)
    mean = _leaf_masses(pt)[leaf] * (2 ** M) * g0
    if ndraws > 0:
        D = np.empty((ndraws, len(xgrid)))
        for d in range(ndraws):
            D[d] = _leaf_masses(pt, draw=True, rng=rng)[leaf] * (2 ** M) * g0
        return mean, D
    return mean, None


# --------------------------------------------------------------------------- #
#  goodness of fit: Bayes factor against the centring family G0                 #
# --------------------------------------------------------------------------- #

def polyatree_logBF(pt):
    """log Bayes factor for H1 (Polya tree) vs H0 (data ~ G0 exactly), summed over all nodes.
    log BF > 0 favours a departure from G0; > ~2.3 (BF>10) is strong evidence against it."""
    M, c = pt["M"], pt["c"]; lbf = 0.0
    for m in range(M):
        a = c * (m + 1) ** 2; nc = pt["counts"][m + 1]
        n0 = nc[0::2]; n1 = nc[1::2]; n = n0 + n1
        lbf += np.sum(n * np.log(2) + betaln(a + n0, a + n1) - betaln(a, a))
    return lbf


# --------------------------------------------------------------------------- #
#  simulation helpers                                                           #
# --------------------------------------------------------------------------- #

def simulate_mixture(n, weights, means, sds, rng):
    weights = np.asarray(weights, float); weights /= weights.sum()
    k = rng.choice(len(weights), size=n, p=weights)
    return rng.normal(np.asarray(means)[k], np.asarray(sds)[k])

def mixture_pdf(grid, weights, means, sds):
    weights = np.asarray(weights, float); weights /= weights.sum()
    g = np.asarray(grid, float); out = np.zeros_like(g)
    for w, m, s in zip(weights, means, sds):
        out += w * np.exp(-0.5 * ((g - m) / s) ** 2) / (s * np.sqrt(2 * np.pi))
    return out

References