Polya Trees¶
A prior on a distribution — and a test of whether a parametric family fits¶
The Dirichlet process is a prior over discrete distributions. A Polya tree (Ferguson 1974; Lavine 1992) is a prior over continuous ones — a random density — built by recursively splitting the sample space and randomising how mass flows down each split.
Map the data to $[0,1]$ through the CDF of a centring distribution $G_0$, $u_i=G_0(x_i)$, and split $[0,1]$ dyadically. At every split the fraction of the parent's mass going left is a Beta variable $$Y_\varepsilon\sim\text{Beta}(\alpha_m,\alpha_m),\qquad \alpha_m=c\,m^2,$$ independent across splits. The choice $\alpha_m=c\,m^2$ makes the random measure a genuine density; the concentration $c$ says how tightly it hugs $G_0$ ($c\to\infty$ returns $G_0$; small $c$ allows wild departures). Since the splits sit at $G_0$'s quantiles, the prior is centred on $G_0$. And it is conjugate — with child counts $n_{\varepsilon0},n_{\varepsilon1}$, $$Y_\varepsilon\mid\text{data}\sim\text{Beta}(\alpha_m+n_{\varepsilon0},\,\alpha_m+n_{\varepsilon1}),$$ so the posterior density is closed-form: no MCMC, just Beta draws for the bands.
That conjugacy also gives a clean goodness-of-fit test: centre on a fitted parametric $G_0$, and if the family is right the data are uniform on the $u$-scale (every $Y\approx\tfrac12$); if not, the branch counts pull the $Y$'s away and a closed-form Bayes factor flags it. We estimate densities (against a kernel estimate), test normality on simulated data (against Kolmogorov–Smirnov, Anderson–Darling and Shapiro–Wilk), and finish on daily stock returns — the non-normality that runs through the whole asset-risk arc.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy import stats
from scipy.stats import gaussian_kde
import polyatree as PT
rng = np.random.default_rng(1)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; PURP="#6b46c1"; GREY="#718096"
print("A Polya tree is a prior over continuous distributions -- and, conjugate, a test of a parametric family.")
A Polya tree is a prior over continuous distributions -- and, conjugate, a test of a parametric family.
1. What the prior looks like¶
Before any data, a Polya tree centred on $N(0,1)$ is a random density scattered around the standard normal. The concentration $c$ controls the scatter: large $c$ keeps draws close to the normal, small $c$ lets them wander. Sampling the prior (all branch counts zero) shows this.
G=stats.norm(0,1); grid=np.linspace(-3.5,3.5,300)
fig,ax=plt.subplots(1,3,figsize=(13,3.8))
for a,cc in zip(ax,[0.5,3.0,30.0]):
pt=dict(counts=[np.zeros(2**m) for m in range(9)], c=cc, M=8, cdf=G.cdf, pdf=G.pdf)
for _ in range(6):
d1=PT.polyatree_density(pt, grid, ndraws=1, rng=rng)[1][0]
a.plot(grid, d1, color=BLUE, lw=.9, alpha=.6)
a.plot(grid, G.pdf(grid), color="k", lw=2, ls="--")
a.set_title(f"prior draws, c={cc:g}"); a.set_xlabel("x"); a.set_ylim(0,0.75)
plt.tight_layout(); plt.show()
print("Each blue curve is a density drawn from the prior; the dashed line is the centring N(0,1). Small c gives")
print("rough, far-ranging densities; large c pins them to the normal. c is the 'trust' placed in the parametric centre.")
Each blue curve is a density drawn from the prior; the dashed line is the centring N(0,1). Small c gives rough, far-ranging densities; large c pins them to the normal. c is the 'trust' placed in the parametric centre.
2. Density estimation — catching a shape the centre misses¶
Data from a bimodal mixture, with the Polya tree centred on a single fitted normal. The posterior density should pull away from the normal and recover the two modes — and match a frequentist kernel density estimate, while additionally supplying a credible band.
w,mu,sd=[0.5,0.5],[-2.2,2.2],[0.7,0.7]
x=PT.simulate_mixture(700,w,mu,sd,rng); m,s=x.mean(),x.std(); G0=stats.norm(m,s)
pt=PT.polyatree_fit(x,G0.cdf,G0.pdf,c=1.0,M=8)
grid=np.linspace(-5.5,5.5,400); dens,D=PT.polyatree_density(pt,grid,ndraws=400,rng=rng)
lo,hi=np.percentile(D,[2.5,97.5],axis=0)
kde=gaussian_kde(x)(grid); truep=PT.mixture_pdf(grid,w,mu,sd)
fig,ax=plt.subplots(figsize=(9,4.4))
ax.hist(x,bins=40,density=True,color=GREY,alpha=.3)
ax.fill_between(grid,lo,hi,color=BLUE,alpha=.2,label="Polya tree 95%")
ax.plot(grid,dens,color=BLUE,lw=2.3,label="Polya tree density")
ax.plot(grid,kde,color=GREEN,lw=1.7,ls="-.",label="kernel density (frequentist)")
ax.plot(grid,G0.pdf(grid),color=RED,lw=1.7,ls="--",label="centring normal G0")
ax.plot(grid,truep,color="k",lw=1.2,ls=":",label="true density")
ax.set_xlabel("x"); ax.set_ylabel("density"); ax.set_title("Polya tree recovers the bimodality the normal centre misses"); ax.legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
print(f"L1 error vs truth: Polya tree {np.trapezoid(np.abs(dens-truep),grid):.3f}, KDE {np.trapezoid(np.abs(kde-truep),grid):.3f}.")
print("The tree starts from the normal but the data drag it to two modes -- it and the kernel estimate agree, and")
print("the tree adds a credible band. The centring normal (red) is exactly the null the next section tests.")
L1 error vs truth: Polya tree 0.191, KDE 0.309. The tree starts from the normal but the data drag it to two modes -- it and the kernel estimate agree, and the tree adds a credible band. The centring normal (red) is exactly the null the next section tests.
3. Goodness of fit — is it normal?¶
Centre the tree on the fitted normal and read the Bayes factor against it: $\log\text{BF}>0$ is evidence of a departure, $>\!2.3$ (BF > 10) strong. We run it on data that is normal (should pass) and on right-skewed data (should reject), and place it beside the classical normality tests — Kolmogorov–Smirnov, Anderson–Darling, Shapiro–Wilk — which return $p$-values for the same question.
def normality_report(x, label):
m,s=x.mean(),x.std(); G0=stats.norm(m,s)
pt=PT.polyatree_fit(x,G0.cdf,G0.pdf,c=1.0,M=8); lbf=PT.polyatree_logBF(pt)
import warnings
with warnings.catch_warnings():
# scipy.stats.anderson warns that a p-value `method` will become required in 1.17.
# We use only the STATISTIC, never a p-value, so the change does not affect this table.
warnings.simplefilter("ignore", FutureWarning)
ks=stats.kstest(x,G0.cdf).pvalue; sh=stats.shapiro(x).pvalue; ad=stats.anderson(x,"norm").statistic
verdict="REJECT normal" if lbf>2.3 else "normal OK"
print(f"{label:16s} logBF {lbf:8.1f} -> {verdict:14s} | KS p {ks:.1e} Shapiro p {sh:.1e} AD stat {ad:.1f}")
return pt, G0
print("data Polya-tree Bayes factor | frequentist normality tests")
ptn,Gn=normality_report(rng.standard_normal(1000), "truly normal")
pts,Gs=normality_report(rng.gamma(2.0,1.0,1000), "skewed (gamma)")
ptt,Gt=normality_report(stats.t(df=3).rvs(1000,random_state=rng), "heavy-tailed t3")
grid=np.linspace(-6,10,400)
fig,ax=plt.subplots(1,2,figsize=(12.5,4.2))
for a,(pt,G0,xx,ttl) in zip(ax,[(pts,Gs,None,"skewed data"),(ptt,Gt,None,"heavy-tailed data")]):
dn,_=PT.polyatree_density(pt,grid); a.plot(grid,dn,color=BLUE,lw=2.2,label="Polya tree"); a.plot(grid,G0.pdf(grid),color=RED,lw=1.7,ls="--",label="fitted normal")
a.set_title(ttl); a.set_xlabel("x"); a.set_ylabel("density"); a.legend(frameon=False,fontsize=8)
ax[0].set_xlim(-3,10); ax[1].set_xlim(-6,6)
plt.tight_layout(); plt.show()
print("The Bayes factor and the classical p-values agree: normal data passes (logBF<0), skew and heavy tails are")
print("rejected decisively. But the Polya tree does more than reject -- it hands back the density that fits instead.")
data Polya-tree Bayes factor | frequentist normality tests truly normal logBF -15.5 -> normal OK | KS p 9.3e-01 Shapiro p 5.6e-01 AD stat 0.2 skewed (gamma) logBF 88.8 -> REJECT normal | KS p 5.0e-08 Shapiro p 2.4e-23 AD stat 19.3 heavy-tailed t3 logBF 28.9 -> REJECT normal | KS p 6.6e-05 Shapiro p 5.5e-21 AD stat 12.8
The Bayes factor and the classical p-values agree: normal data passes (logBF<0), skew and heavy tails are rejected decisively. But the Polya tree does more than reject -- it hands back the density that fits instead.
4. Are daily stock returns normal?¶
The asset-risk arc leans on the answer being no — GARCH volatility, Student-$t$ innovations, extreme-value tails and copulas all exist because returns are not Gaussian. The Polya tree tests it directly: centre on the fitted normal for 15 years of daily S&P 500 (SPY) returns and read the Bayes factor and the density.
d=pd.read_csv("spy_returns.csv"); r=d["ret"].to_numpy(); m,s=r.mean(),r.std(); G0=stats.norm(m,s)
pt=PT.polyatree_fit(r,G0.cdf,G0.pdf,c=1.0,M=9); lbf=PT.polyatree_logBF(pt)
print(f"SPY daily returns (n={len(r)}): skew {stats.skew(r):.2f}, excess kurtosis {stats.kurtosis(r):.1f}")
print(f"Polya-tree log Bayes factor vs normal = {lbf:.0f} (BF astronomically large -> normality crushed)")
print(f"frequentist: KS p {stats.kstest(r,G0.cdf).pvalue:.1e}, Shapiro p {stats.shapiro(r[:4000]).pvalue:.1e}")
grid=np.linspace(-6,6,500); dens,D=PT.polyatree_density(pt,grid,ndraws=300,rng=rng); lo,hi=np.percentile(D,[2.5,97.5],axis=0)
fig,ax=plt.subplots(1,2,figsize=(12.5,4.3))
ax[0].hist(r,bins=120,density=True,color=GREY,alpha=.35)
ax[0].fill_between(grid,lo,hi,color=BLUE,alpha=.2); ax[0].plot(grid,dens,color=BLUE,lw=2,label="Polya tree"); ax[0].plot(grid,G0.pdf(grid),color=RED,lw=1.8,ls="--",label="fitted normal")
ax[0].set_xlim(-6,6); ax[0].set_xlabel("daily return (%)"); ax[0].set_ylabel("density"); ax[0].set_title("SPY return density: Polya tree vs normal"); ax[0].legend(frameon=False)
ax[1].fill_between(grid,lo,hi,color=BLUE,alpha=.2); ax[1].plot(grid,dens,color=BLUE,lw=2); ax[1].plot(grid,G0.pdf(grid),color=RED,lw=1.8,ls="--")
ax[1].set_yscale("log"); ax[1].set_xlim(-6,6); ax[1].set_ylim(1e-4,1); ax[1].set_xlabel("daily return (%)"); ax[1].set_ylabel("density (log)"); ax[1].set_title("Log scale: the fat tails and sharp peak")
plt.tight_layout(); plt.show()
vals, cnts = np.unique(r, return_counts=True); n_rep = int((cnts - 1).sum())
print("\nOne assumption worth checking rather than assuming, since this is a goodness-of-fit notebook:")
print("KS and Shapiro both assume a CONTINUOUS distribution, i.e. no repeated values. This series has")
print("%d repeated observations out of %d -- all of them exact zeros, days the index closed unchanged." % (n_rep, len(r)))
print("(R's ks.test warns about exactly this; scipy stays silent.) At p = %.0e the verdict is nowhere near" % stats.kstest(r, G0.cdf).pvalue)
print("sensitive to %d ties, so nothing here changes -- but note which method had to make the assumption." % n_rep)
print("The Polya tree partitions the line and counts what falls in each cell, so repeated values violate")
print("nothing it relies on. That is a small structural advantage, and it is the kind that matters on")
print("price data, where exact zeros and rounded ticks are ordinary rather than pathological.\n")
print("The Polya-tree density has a taller, sharper peak and much fatter tails than the normal -- the log-scale panel")
print("shows the tails sitting orders of magnitude above the Gaussian. This is exactly the leptokurtosis the asset-risk")
print("arc models with Student-t GARCH, EVT and copulas; here a nonparametric prior confirms and quantifies it directly.")
SPY daily returns (n=3772): skew -0.72, excess kurtosis 11.5 Polya-tree log Bayes factor vs normal = 257 (BF astronomically large -> normality crushed) frequentist: KS p 5.1e-33, Shapiro p 3.7e-44
One assumption worth checking rather than assuming, since this is a goodness-of-fit notebook: KS and Shapiro both assume a CONTINUOUS distribution, i.e. no repeated values. This series has 11 repeated observations out of 3772 -- all of them exact zeros, days the index closed unchanged. (R's ks.test warns about exactly this; scipy stays silent.) At p = 5e-33 the verdict is nowhere near sensitive to 11 ties, so nothing here changes -- but note which method had to make the assumption. The Polya tree partitions the line and counts what falls in each cell, so repeated values violate nothing it relies on. That is a small structural advantage, and it is the kind that matters on price data, where exact zeros and rounded ticks are ordinary rather than pathological. The Polya-tree density has a taller, sharper peak and much fatter tails than the normal -- the log-scale panel shows the tails sitting orders of magnitude above the Gaussian. This is exactly the leptokurtosis the asset-risk arc models with Student-t GARCH, EVT and copulas; here a nonparametric prior confirms and quantifies it directly.
5. Summary¶
A Polya tree is a prior over continuous distributions: split the space dyadically, randomise the mass flow with Beta variables, and centre the whole thing on a base $G_0$. Because it is conjugate, the posterior density is closed-form — the only model in this arc that needs no sampler at all, just Beta draws for the bands. Centred on a fitted parametric family it becomes a goodness-of-fit test: the Bayes factor against $G_0$ passed truly normal data, rejected skew and heavy tails, and — like the classical Kolmogorov–Smirnov / Anderson–Darling / Shapiro–Wilk tests — crushed the normality of daily S&P 500 returns, while additionally returning the fat-tailed density that fits instead.
That last point ties the whole portfolio together: the asset-risk arc's Student-$t$ GARCH, extreme-value tails and copulas all exist because returns are non-normal, and here a Bayesian nonparametric prior confirms it directly and quantifies the leptokurtosis. Unlike the Dirichlet process (discrete distributions), the Polya tree is a prior on a density, and its conjugacy makes it the natural Bayesian companion to the frequentist distribution tests.
This closes the Bayesian nonparametrics arc — Dirichlet-process mixtures, the hierarchical DP, Gaussian-process regression and classification, penalised splines, nonparametric hazards, and now Polya trees — priors placed, in turn, on clusters, functions, hazards, and whole distributions.