Gaussian vs Cauchy model, common prior

Closed-form vs SIR-based inference

1 Models

Let us assume that observations \(x_1,\ldots,x_n\) are independent and identifically distributed and that we are interested in comparing competing models

1. Model \({\cal M}_1: \ Gaussian(\mu,1)\), or
2. Model \({\cal M}_2: \ Cauchy(\mu,1)\),

where \(\mu\) is the mean of the \(X\)’s under the Gaussian model, \({\cal M}_1\), but not under the Cauchy model, \({\cal M}_2\), since the Cauchy distribution has undefined mean. Nonetheless, \(\mu\) can still be interpreted as the median/model of the \(X\)’s under the Cauchy model.

2 Prior

Since \(\mu\) represent the same quantity (location) under both models, we will assume a common prior for \(\mu\): \[ \mu \sim N(0,\tau_0^2), \] where, for simplicity, we will assume that \(\tau_0^2=10\).

3 Posterior

Under the Gaussian model, the posterior distribution of \(\mu\) follows from simple conjugacy: \[ \mu|x_{1:n},{\cal M}_1 \sim N(\mu_1,\tau_1^2), \] where \(\mu_1=\sum_{i=1}^n x_i/(n+1/\tau_0^2)\) and \(1/(n+1/\tau_0^2)\). So, posterior inference can be performed in closed form.

However, this is not the case for the posterior of \(\mu\) under the Cauchy model: \[ p(\mu|x_{1:n},{\cal M}_1) \propto \left[\prod_{i=1}^n p(x_i|\mu,{\cal M}_2)\right] p(\mu) \propto \left[\prod_{i=1}^n \frac{1}{1+(x_i-\mu)^2}\right]\exp\{-0.5\mu^2/\tau_0^2\}, \] which does not follow any known family of distributions.

# Normalized posterior for model 1
post1 = function(sumx,mu,tau2){
  dnorm(mu,sumx/(n+1/tau2),sqrt(1/(n+1/tau2)))
}
# Unnormalized posterior for model 2
post2 = function(x,mu,tau2){
  prod(dt(mu-x,df=1))*dnorm(mu,0,sqrt(tau2))
}
# Unnormalized posterior for model 2 (log)
logpost2 = function(x,mu,tau2){
  sum(dt(mu-x,df=1,log=TRUE))+dnorm(mu,0,sqrt(tau2),log=TRUE)
}

4 Turning the Bayesian crank

# Hyperparameters
tau2=10

# Observed data
x    = c(-4.3,3.2)
sumx = sum(x)
n    = length(x)

# Plotting the posterior densities
par(mfrow=c(1,2))

mu = seq(-5,5,length=1000)
plot(mu,post1(sumx,mu,tau2),type="l",lwd=2,ylab="Posterior")
title("Model 1: Gaussian")
for (i in 1:n)
  points(x[i],0,pch=16,col=3)

mu = seq(-10,10,length=1000)
lp = rep(0,1000)
for (i in 1:1000)
  lp[i] = logpost2(x,mu[i],tau2)
lp1 = exp(lp-max(lp))
plot(mu,lp1,type="l",lwd=2,ylab="Prior*likelihood")
title("Model 2: Cauchy")
for (i in 1:n)
  points(x[i],0,pch=16,col=3)

5 Computing \(p(x_{1:n}|{\cal M}_2)\) via simple Reimann integration

M      = 1000000
mus   = seq(-15,15,length=M)
h       = mus[2]-mus[1]
prod = rep(0,M)
for (i in 1:M)
prod[i] = h*post2(x,mus[i],tau2)
kappa = sum(prod)
kappa

## [1] 0.000775713

We can now plot both posterior densities under the same figure:

par(mfrow=c(1,1))
plot(mu,post1(sumx,mu,tau2),type="l",lwd=2,xlab=expression(mu),
ylab="Posterior density")
lines(mu,exp(lp-log(kappa)),col=2,lwd=2)
lines(mu,dnorm(mu,0,sqrt(10)),col=4,lwd=2)
legend("topright",legend=c("Prior of mu","Posterior of mu (Gaussian model)","Posterior of mu (Cauchy model)"),col=c(4,1:2),lwd=2,lty=1)
for (i in 1:n)
  points(x[i],0,pch=16,col=3)

6 Computing other posterior summaries

Below we compute, via simple Reimann integration, the following summaries:

1. \(E(\mu|x_{1:n},{\cal M}_2)\)
2. \(E(\mu^2|x_{1:n},{\cal M}_2)\)
3. \(Pr(\mu<1|x_{1:n},{\cal M}_2)\)

M      = 1000000
mus    = seq(-15,15,length=M)
h      = mus[2]-mus[1]
prod   = matrix(0,M,3)
for (i in 1:M){
  prod[i,1] = mus[i]*h*post2(x,mus[i],tau2)/kappa
  prod[i,2] = mus[i]^2*h*post2(x,mus[i],tau2)/kappa
  prod[i,3] = ifelse(mus[i]<1,1,0)*h*post2(x,mus[i],tau2)/kappa
}
E = sum(prod[,1])
E2 = sum(prod[,2])
V = E2-E^2
SD = sqrt(V)

tab = round(rbind(c(sumx/(n+1/tau2),sqrt(1/(n+1/tau2)),
pnorm(1,sumx/(n+1/tau2),sqrt(1/(n+1/tau2)))),
c(E,SD,sum(prod[,3]))),3)
rownames(tab) = c("Gaussian model","Cauchy model")
colnames(tab) = c("E(mu)","SD(mu)","Pr(mu<1)")
tab

##                 E(mu) SD(mu) Pr(mu<1)
## Gaussian model -0.524   0.69    0.986
## Cauchy model   -0.066   2.95    0.550

7 SIR-based posterior inference under the Cauchy model

N = 1000000
mus = rnorm(N,0,sqrt(10))
w = rep(0,N)
for (i in 1:N)
  w[i] = prod(dt(x-mus[i],df=1))
ind = sample(1:N,size=N,prob=w,replace=TRUE)
mus1 = mus[ind]

hist(mus1,prob=TRUE,xlab=expression(mu),ylab="Posterior density",main="",ylim=c(0,0.25))
mu = seq(-10,10,length=1000)
lines(mu,exp(lp-log(kappa)),col=2,lwd=2)
legend("topleft",legend=c("Exact posterior","MC-based posterior draws"),col=2:1,lty=1,lwd=2)

7.1 Summarizing the results

tab = rbind(tab,round(c(mean(mus1),sqrt(var(mus1)),mean(mus1<1)),3))
rownames(tab) = c("Gaussian model (closed form)","Cauchy model (Reimann)","Cauchy model (SIR)")
tab = cbind(tab[c(1,3),],
rbind(c(qnorm(0.05,sumx/(n+1/tau2),sqrt(1/(n+1/tau2))),
qnorm(0.95,sumx/(n+1/tau2),sqrt(1/(n+1/tau2)))),
quantile(mus1,c(0.05,0.95))))
tab

##                               E(mu) SD(mu) Pr(mu<1)        5%       95%
## Gaussian model (closed form) -0.524  0.690    0.986 -1.658866 0.6112473
## Cauchy model (SIR)           -0.065  2.952    0.550 -4.635423 3.8669082

8 Model comparison

The prior predictive for the Gaussian model is easy to derive \[ p(x_{1:n}|{\cal M}_1) = \prod_{i=1}^n p_N(x_i|0,1+\tau^2). \]

priorpred.M1 =prod(dnorm(x,0,tau2+1))
priorpred.M2.Reimann = kappa
priorpred.M2.MonteCarlo = mean(w)
BF.Reimann = priorpred.M1/priorpred.M2.Reimann
BF.MonteCarlo = priorpred.M1/priorpred.M2.MonteCarlo
c(priorpred.M1,priorpred.M2.Reimann,priorpred.M2.MonteCarlo)

## [1] 0.0011680885 0.0007757130 0.0007756855

c(BF.Reimann,BF.MonteCarlo)

## [1] 1.505826 1.505879