Polya Trees (R)¶
From-scratch nonparametric density prior, with the classical normality tests¶
The R counterpart to polyatree_python.ipynb. A Polya tree centres a prior on a distribution $G_0$: map data to $[0,1]$ by $u=G_0(x)$, split dyadically, and randomise each split's mass with $Y_\varepsilon\sim\text{Beta}(cm^2,cm^2)$. It is conjugate — $Y_\varepsilon\mid\text{data}\sim\text{Beta}(cm^2+n_{\varepsilon0},\,cm^2+n_{\varepsilon1})$ — so the posterior density is closed-form (no MCMC). We port it to base R, estimate a density, and use it as a goodness-of-fit test against the classical Kolmogorov–Smirnov (ks.test) and Shapiro–Wilk (shapiro.test) normality tests, finishing on daily S&P 500 returns.
options(repr.plot.width=12, repr.plot.height=4.2)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
BLUE<-"#2b6cb0"; RED<-"#c53030"; GREEN<-"#2f855a"; GREY<-"#718096"
pt_fit<-function(x, cdf, pdf, c=1, M=8){ u<-pmin(pmax(cdf(x),1e-12),1-1e-12)
counts<-lapply(0:M, function(m){ idx<-pmin(floor(u*2^m),2^m-1); tabulate(idx+1, nbins=2^m) })
list(counts=counts, c=c, M=M, cdf=cdf, pdf=pdf) }
leaf_masses<-function(pt, draw=FALSE){ M<-pt$M; c<-pt$c; mass<-1
for(m in 0:(M-1)){ a<-c*(m+1)^2; nc<-pt$counts[[m+2]]; n0<-nc[seq(1,length(nc),2)]; n1<-nc[seq(2,length(nc),2)]
y<-if(draw) rbeta(length(n0),a+n0,a+n1) else (a+n0)/(2*a+n0+n1)
nm<-numeric(2^(m+1)); nm[seq(1,2^(m+1),2)]<-mass*y; nm[seq(2,2^(m+1),2)]<-mass*(1-y); mass<-nm }
mass }
pt_density<-function(pt, xg, ndraws=0){ M<-pt$M; u<-pmin(pmax(pt$cdf(xg),1e-12),1-1e-12); leaf<-pmin(floor(u*2^M),2^M-1)+1; g0<-pt$pdf(xg)
meand<-leaf_masses(pt)[leaf]*2^M*g0
if(ndraws>0){ D<-matrix(0,ndraws,length(xg)); for(d in 1:ndraws) D[d,]<-leaf_masses(pt,TRUE)[leaf]*2^M*g0; return(list(mean=meand,draws=D)) }
list(mean=meand,draws=NULL) }
pt_logBF<-function(pt){ M<-pt$M; c<-pt$c; lbf<-0
for(m in 0:(M-1)){ a<-c*(m+1)^2; nc<-pt$counts[[m+2]]; n0<-nc[seq(1,length(nc),2)]; n1<-nc[seq(2,length(nc),2)]; n<-n0+n1
lbf<-lbf+sum(n*log(2)+lbeta(a+n0,a+n1)-lbeta(a,a)) }
lbf }
cat("from-scratch base-R Polya tree (density + goodness-of-fit Bayes factor) ready\n")
from-scratch base-R Polya tree (density + goodness-of-fit Bayes factor) ready
1. The prior — random densities around a normal centre¶
Before data, a Polya tree centred on $N(0,1)$ is a random density scattered around the standard normal; the concentration $c$ sets the scatter (large $c$ hugs the normal, small $c$ wanders).
grid<-seq(-3.5,3.5,length=300)
par(mfrow=c(1,3), mar=c(4,4,3,1))
for(cc in c(0.5,3,30)){ pt<-list(counts=lapply(0:8,function(m) numeric(2^m)), c=cc, M=8, cdf=pnorm, pdf=dnorm)
plot(grid, dnorm(grid), type="n", ylim=c(0,0.75), xlab="x", ylab="density", main=bquote(c==.(cc)))
for(i in 1:6){ set.seed(i*10+cc); lines(grid, pt_density(pt,grid,ndraws=1)$draws[1,], col=adjustcolor(BLUE,.6), lwd=.9) }
lines(grid, dnorm(grid), col="black", lwd=2, lty=2) }
par(mfrow=c(1,1))
cat("Each blue curve is a prior density draw; dashed = the centring N(0,1). Small c wanders, large c pins to normal.\n")
Each blue curve is a prior density draw; dashed = the centring N(0,1). Small c wanders, large c pins to normal.
2. Density estimation vs a kernel estimate¶
Bimodal data, tree centred on a fitted normal: the posterior density pulls away from the normal to the two modes and matches base R's kernel density(), with a credible band.
set.seed(1); k<-sample(1:2,700,replace=TRUE,prob=c(.5,.5)); x<-rnorm(700, c(-2.2,2.2)[k], 0.7)
m<-mean(x); s<-sd(x); G0cdf<-function(t) pnorm(t,m,s); G0pdf<-function(t) dnorm(t,m,s)
pt<-pt_fit(x,G0cdf,G0pdf,c=1,M=8); grid<-seq(-5.5,5.5,length=400); pd<-pt_density(pt,grid,ndraws=400)
lo<-apply(pd$draws,2,quantile,.025); hi<-apply(pd$draws,2,quantile,.975); kd<-density(x,n=400,from=-5.5,to=5.5)
par(mar=c(4,4,3,1)); hist(x,breaks=40,freq=FALSE,col=adjustcolor(GREY,.3),border=NA,xlim=c(-5.5,5.5),ylim=c(0,max(hi)),xlab="x",main="Polya tree recovers the bimodality the normal centre misses")
polygon(c(grid,rev(grid)),c(lo,rev(hi)),col=adjustcolor(BLUE,.2),border=NA); lines(grid,pd$mean,col=BLUE,lwd=2.4)
lines(kd,col=GREEN,lwd=1.7,lty=4); lines(grid,G0pdf(grid),col=RED,lwd=1.7,lty=2)
legend("topright",c("Polya tree","kernel density()","centring normal"),col=c(BLUE,GREEN,RED),lwd=2,lty=c(1,4,2),bty="n",cex=.8)
cat("The tree starts from the normal but the data drag it to two modes, agreeing with the kernel estimate; the\n")
cat("centring normal (red) is the null the goodness-of-fit test rejects next.\n")
The tree starts from the normal but the data drag it to two modes, agreeing with the kernel estimate; the
centring normal (red) is the null the goodness-of-fit test rejects next.
3. Goodness of fit — Bayes factor vs ks.test / shapiro.test¶
Centre on the fitted normal and read the log Bayes factor ($>2.3$ = strong rejection), beside the classical normality $p$-values, on normal / skewed / heavy-tailed data.
report<-function(x,label){ m<-mean(x); s<-sd(x); pt<-pt_fit(x,function(t)pnorm(t,m,s),function(t)dnorm(t,m,s),c=1,M=8)
lbf<-pt_logBF(pt); ks<-ks.test(x,"pnorm",m,s)$p.value; sh<-shapiro.test(x)$p.value
cat(sprintf("%-16s logBF %8.1f -> %-14s | KS p %.1e Shapiro p %.1e\n", label, lbf, ifelse(lbf>2.3,"REJECT normal","normal OK"), ks, sh))
pt }
set.seed(2)
cat("data Polya-tree Bayes factor | frequentist tests\n")
p1<-report(rnorm(1000), "truly normal")
p2<-report(rgamma(1000,2,1), "skewed (gamma)")
p3<-report(rt(1000,df=3), "heavy-tailed t3")
cat("\nThe Bayes factor and the classical p-values agree: normal passes, skew and heavy tails are rejected. The tree\n")
cat("adds what the tests cannot -- the density that fits instead (shown for the returns below).\n")
data Polya-tree Bayes factor | frequentist tests
truly normal logBF -21.2 -> normal OK | KS p 9.4e-01 Shapiro p 5.2e-01
skewed (gamma) logBF 91.2 -> REJECT normal | KS p 2.2e-09 Shapiro p 1.0e-24
heavy-tailed t3 logBF 27.9 -> REJECT normal | KS p 3.4e-05 Shapiro p 7.0e-22
The Bayes factor and the classical p-values agree: normal passes, skew and heavy tails are rejected. The tree
adds what the tests cannot -- the density that fits instead (shown for the returns below).
4. Are daily S&P 500 returns normal?¶
The asset-risk arc assumes not — Student-$t$ GARCH, extreme-value tails and copulas all address non-normal returns. Test it directly on 15 years of daily SPY returns.
d<-read.csv("spy_returns.csv"); r<-d$ret; m<-mean(r); s<-sd(r)
pt<-pt_fit(r,function(t)pnorm(t,m,s),function(t)dnorm(t,m,s),c=1,M=9); lbf<-pt_logBF(pt)
cat(sprintf("SPY daily returns (n=%d): excess kurtosis %.1f\n", length(r), mean((r-m)^4)/s^4-3))
cat(sprintf("Polya-tree log Bayes factor vs normal = %.0f (normality crushed); KS p %.1e, Shapiro p %.1e\n",
lbf, ks.test(r,"pnorm",m,s)$p.value, shapiro.test(r[1:4000])$p.value))
grid<-seq(-6,6,length=500); pd<-pt_density(pt,grid,ndraws=250); lo<-apply(pd$draws,2,quantile,.025); hi<-apply(pd$draws,2,quantile,.975)
par(mfrow=c(1,2), mar=c(4,4,3,1))
hist(r,breaks=120,freq=FALSE,col=adjustcolor(GREY,.35),border=NA,xlim=c(-6,6),xlab="daily return (%)",main="SPY density: Polya tree vs normal")
polygon(c(grid,rev(grid)),c(lo,rev(hi)),col=adjustcolor(BLUE,.2),border=NA); lines(grid,pd$mean,col=BLUE,lwd=2); lines(grid,dnorm(grid,m,s),col=RED,lwd=1.8,lty=2)
legend("topright",c("Polya tree","fitted normal"),col=c(BLUE,RED),lwd=2,lty=c(1,2),bty="n")
plot(grid,pd$mean,type="l",col=BLUE,lwd=2,log="y",ylim=c(1e-4,1),xlim=c(-6,6),xlab="daily return (%)",ylab="density (log)",main="Log scale: fat tails and sharp peak")
lines(grid,dnorm(grid,m,s),col=RED,lwd=1.8,lty=2)
par(mfrow=c(1,1))
n_rep <- sum(duplicated(r))
cat(sprintf("\nThe ks.test warning above is not noise: KS and Shapiro assume a CONTINUOUS distribution, and\n"))
cat(sprintf("this series has %d repeated values out of %d -- all exact zeros, days the index closed unchanged.\n",
n_rep, length(r)))
cat("At these p-values the conclusion is nowhere near sensitive to that, but note which method needed the\n")
cat("assumption. The Polya tree partitions the line and counts what lands in each cell, so repeated values\n")
cat("violate nothing it relies on -- a small structural advantage, and one that matters on price data where\n")
cat("exact zeros and rounded ticks are ordinary rather than pathological.\n")
cat("\nTaller sharper peak, much fatter tails: on the log scale the Polya-tree tails sit orders of magnitude above the\n")
cat("normal. The nonparametric prior confirms and quantifies the leptokurtosis the asset-risk arc models throughout.\n")
SPY daily returns (n=3772): excess kurtosis 11.5
Warning message in ks.test.default(r, "pnorm", m, s): "ties should not be present for the one-sample Kolmogorov-Smirnov test"
Polya-tree log Bayes factor vs normal = 257 (normality crushed); KS p 6.2e-33, Shapiro p 3.7e-44
The ks.test warning above is not noise: KS and Shapiro assume a CONTINUOUS distribution, and
this series has 11 repeated values out of 3772 -- all exact zeros, days the index closed unchanged.
At these p-values the conclusion is nowhere near sensitive to that, but note which method needed the
assumption. The Polya 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.
Taller sharper peak, much fatter tails: on the log scale the Polya-tree tails sit orders of magnitude above the
normal. The nonparametric prior confirms and quantifies the leptokurtosis the asset-risk arc models throughout.
5. Summary¶
The from-scratch base-R Polya tree reproduces the Python results: a conjugate prior over continuous distributions whose posterior density is closed-form, matching base R's kernel density() on bimodal data and — centred on a fitted normal — serving as a goodness-of-fit test. Its Bayes factor agreed with ks.test and shapiro.test: normal data passed, skew and heavy tails were rejected, and daily S&P 500 returns were decisively non-normal, with a density showing the sharp peak and fat tails the arc's Student-$t$ GARCH, EVT and copula models exist to capture.
Together with polyatree_python.ipynb (the prior-scatter picture, the density-vs-KDE recovery, and the full test battery) this is the Polya tree — a Bayesian nonparametric prior placed directly on a distribution, and the natural companion to the classical distribution tests. It closes the Bayesian-nonparametrics arc: priors on clusters (Dirichlet process), on functions (Gaussian processes), on hazards (Gamma process), and now on whole densities.