Sequential Bayesian Learning

Beta-Binomial and Pareto-Uniform

1 Beta-Binomial model

The data is a sequence of Bernoulli trials, i.e. \(x_1,\ldots,x_n\) are iid \(Bernoulli(\theta)\) for \(\theta \in (0,1)\). The prior for \(\theta\) is uniform, i.e. \(p(\theta)=1\), for \(\theta \in (0,1)\).

Then, the posterior can be derived as \[ p(\theta|x_{1:n}) \propto \prod_{i=1}^n \theta^{x_i}(1-\theta)^{1-x_i} = \theta^{s_n}(1-\theta)^{n-s_n}, \] where \(x_{1:n}=\{x_1,\ldots,x_n\}\) and \(s_n=\sum_{i=1}^n x_i\). This posterior density exhibits the kernel of a Beta distribution with parameters \(\alpha_n=1+s_n\) and \(\beta_n=1+n-s_n\), and the following summaries: \[\begin{eqnarray*} E(\theta|x_{1:n}) &=& \frac{1+s_n}{1+n}\\ V(\theta|x_{1:n}) &=& \frac{(1+s_n)(1+n-s_n)}{(2+n)^2(3+n)}\\ Mode(\theta|x_{1:n})&=& \frac{s_n}{n} = {\widehat \theta}_{mle} \end{eqnarray*}\]

Recall that, when \(\theta \sim Beta(a,b)\), then \[\begin{eqnarray*} p(\theta) &\propto& \theta^{a-1}(1-\theta)^{b-1},\\ E(\theta) &=& \frac{a}{a+b},\\ Mode(\theta)&=& \frac{a-1}{a+b-2}, \ \mbox{for} \ a,b>1,\\ V(\theta)&=&\frac{ab}{(a+b)^2(a+b+1)}. \end{eqnarray*}\] See, for instance, the Wikipedia page for more details about the Beta distribution.

1.1 Simulating Bernoulli\((\theta)\) trials

set.seed(1235)
theta = 0.2
n     = 200
x     = rbinom(n,1,theta)
par(mfrow=c(1,1))
plot(x,xlab="Trials",ylab="Outcome")

1.2 Sequential MLE of \(\theta\)

theta.mle = cumsum(x)/(1:n)

par(mfrow=c(1,1))
plot(theta.mle,xlab="Trials",ylab="MLE of theta")
abline(h=theta,col=2)

1.3 Sequential posterior inference

alphan = cumsum(x)+1
betan  = 2:(n+1)-cumsum(x)
E      = alphan/(alphan+betan)
M      = (alphan-1)/(alphan+betan-2)

# Posterior percentiles: 2.5%, 50% and 97.5%
IC = cbind(
  qbeta(0.025,alphan,betan),
  qbeta(0.5,alphan,betan),
  qbeta(0.975,alphan,betan)
)

par(mfrow=c(1,1))
plot(E,ylim=c(0,1),xlab="Sample size",ylab="Posterior quantiles theta",type="b",pch=16,cex=0.5)
lines(IC[,1],type="b",col=2,pch=16,cex=0.5)
lines(IC[,2],type="b",col=2,pch=16,cex=0.5)
lines(IC[,3],type="b",col=2,pch=16,cex=0.5)
points(theta.mle,col=3,pch=16,cex=0.5)
abline(h=theta,col=4,lwd=2)
abline(h=0.1,lty=2)
abline(h=0.3,lty=2)
abline(h=0.5,lty=2)
abline(h=0.7,lty=2)
legend("topright",legend=c("Posterior mean","Posterior mode=MLE","Posterior (2.5,50.9.5) percentiles","True theta"),col=c(1,3,2,4),pch=16,bty="n")

2 Pareto-Uniform model

The data is a sequence of uniform trials, i.e. \(x_1,\ldots,x_n\) are iid \(U[0,\theta]\) for \(\theta>0\). Therefore, \[ p(x_1,\ldots,x_n|\theta) = \prod_{i=1}^n \frac{1}{\theta}I(x_i\leq \theta)= \frac{1}{\theta^n}I(x_{(n)} \leq \theta), \] where \(x_{(n)}=\max\{x_1,\ldots,x_n\}\).

The prior for \(\theta\) is a Pareto distribution with hyperparameters \(\theta_0\) and \(\alpha_0\), denoted by \(\theta \sim Pareto(\theta_0,\alpha_0)\), for \(\theta_0,\alpha_0>0\). Its probability density function is \[ p(\theta|\theta_0,\alpha_0) = \frac{\alpha_0 \theta_0^\alpha}{a\theta^{\alpha_0+1}}, \qquad \mbox{for} \ \theta \geq \theta_0, \] and \(p(\theta|\theta_0,\alpha_0)=0\) for \(\theta<\theta_0\). It can be shown that \[\begin{eqnarray*} E(\theta) &=& \frac{\alpha_0\theta_0}{\alpha_0-1}, \ \alpha_0>1\\ Median(\theta)&=& \theta_0 2^{1/\alpha_0},\\ V(\theta)&=&\frac{\theta_0^2\alpha_0}{(\alpha_0-1)^2(\alpha_0-2)}, \ \alpha_0>2. \end{eqnarray*}\] See, for instance, the Wikipedia page for more details about the Beta distribution.

Consequently, the posterior distribution of \(\theta\) is given by \[ p(\theta|x_{1:n}) \propto \left(\frac{1}{\theta^n}I(\theta \geq x_{(n)})\right) \left(\frac{1}{\theta^{\alpha_0+1}}I(\theta \geq \theta_0)\right) \propto \frac{1}{\theta^{\alpha_1+1}}I(\theta \geq \theta_1), \] where \(\alpha_1=\alpha_0+n\) and \(\theta_1=\max\{x_{(n)},\theta_0\}\). Therefore, the posterior of \(\theta\) is also Pareto, i.e. \(\theta|x_{1:n} \sim Pareto(\theta_1,\alpha_1)\), with \[\begin{eqnarray*} E(\theta|x_{1:n}) &=& \frac{\theta_1\alpha_1}{\alpha_1+1}\\ V(\theta|x_{1:n}) &=& \frac{\theta_1^2 \alpha_1}{(\alpha_1-1)^2(\alpha_1-2)} \end{eqnarray*}\]

2.1 Simulating Uniform\((0,\theta)\) trials

#install.packages("EnvStats")
library("EnvStats")

## 
## Attaching package: 'EnvStats'

## The following objects are masked from 'package:stats':
## 
##     predict, predict.lm

set.seed(165431)
n     = 50
theta = 2
x     = runif(n,0,theta)
plot(x,xlab="Trials",ylab="x")

2.2 Sequential MLE of \(\theta\)

theta.mle = cummax(x)

par(mfrow=c(1,1))
plot(theta.mle,xlab="Trials",ylab="MLE of theta",pch=16,cex=0.5)
abline(h=theta,col=2)
title(paste("max(x)=",round(max(x),3),sep=""))

2.3 Posterior inference

# Prior hyperparameters
theta0 = 1.5
alpha0 = 0.1

# Posterior hyperparameters and posterior summaries
theta1 = max(x,theta0)
alpha1 = alpha0+n
E      = theta1*alpha1/(alpha1-1)
V      = theta1^2*alpha1/((alpha1-1)^2*(alpha1-2))
quant   = qpareto(c(0.05,0.5,0.95),theta1,alpha1)

# Graphical summary
par(mfrow=c(1,1))
thetas = seq(theta1,3,length=1000)
plot(thetas,dpareto(thetas,theta1,alpha1),xlim=range(thetas,x),
     type="l",lwd=2,xlab=expression(theta),ylab="Posterior density")
abline(v=quant[1],lty=2)
abline(v=quant[2],lty=2)
abline(v=quant[3],lty=2)
abline(v=E,col=2)
points(x,rep(0,n),pch=16)
abline(v=theta,col=4)

2.4 Sequential posterior inference

# Sequential learning
stats = matrix(0,n-4,8)
for (i in 5:n){
  xx = x[1:i]
  theta1 = max(xx,theta0)
  alpha1 = alpha0+i
  E = theta1*alpha1/(alpha1-1)
  V = theta1^2*alpha1/((alpha1-1)^2*(alpha1-2))
  quant = qpareto(c(0.1,0.5,0.9),theta1,alpha1)
  stats[i-4,] = c(i,theta1,alpha1,E,V,quant)
}

par(mfrow=c(1,1))
plot(stats[,1],stats[,4],ylim=range(stats[,6:8]),type="l",ylab="Posterior density of theta",xlab="Sample size")
lines(stats[,1],stats[,7],col=2)
lines(stats[,1],stats[,6],col=2,lty=2)
lines(stats[,1],stats[,8],col=2,lty=2)
lines(5:n,cummax(x)[5:n],col=3)
abline(h=theta,col=4)
legend("topright",legend=c("Quantiles 10-50-90","Posterior mean","max(x1,...,xn)"),col=c(2,1,3),lty=1,bty="n")