Data augmentation

Student’s t data

1 Student’s \(t\) data

Let us assume that \(\mbox{data}=\{y_1,\ldots,y_n\}\) are iid \(t_\nu(\mu,\sigma^2)\), where \((\nu,\sigma^2)\) are know. Therefore, our inferential problem is to learn about \(\mu\) conditionally on the data. The prior for \(\mu\) is \(N(\mu_0,\tau_0^2\))$, for known hyperparameters \((\mu_0,\tau_0^2)\).

2 Intractable posterior

It is easy to see that \[ p(\mu|data) \propto \exp\{-0.5(\mu^2+2\mu \mu_0)/\tau_0^2\} \prod_{i=1}^n p_t(y_i|\mu,\sigma^2,\nu), \] is not of known form. Therefore, posterior moments and other summaries will need to be computaed approximately.

3 Augmenting the problem (Gibbs full conditionals)

A standard result in probability is that \(y \sim t_\nu(\mu,\sigma^2)\) can be obtained as a marginal distribution of the following joint distribution \[ p(y,\lambda|\mu,\sigma^2,\nu) = p(y|\lambda,\mu,\sigma^2)p(\lambda|\nu), \] where \(y|\lambda,\mu,\sigma^2 \sim N(\mu,\lambda \sigma^2)\) and \(\lambda | \nu \sim IG(\nu/2,\nu/2)\). In other words, the intractable likelihood becomes tractable once it is augmented by parameters \(\lambda_1,\ldots,\lambda_n\). Therefore, we can sample \(\{(\mu,\lambda_1,\ldots,\lambda_n)^{(j)}\}_{j=1}^M\) iteratively as follows.

Assume that \(\mu^{(0)}={\bar y}_n\), the sample mean. Then, cycle through the full conditionals (Gibbs sampler), for \(j=1,\ldots,M\) iterations:

1. For \(i=1,\ldots,n\), sample \[ \lambda^{(j)}|y_i,\mu,\sigma^2,\nu \sim \left(\frac{\nu+1}{2},\frac{\nu+(y_i-\mu)^2/\sigma^2}{2} \right) \]
2. Sample \(\mu\) from \[ \mu | \mbox{data},\lambda_1,\ldots,\lambda_n,\sigma^2,\nu \sim N(\mu_1,\tau_1^2), \] where \[ \tau_1^{-2} = \tau_0^{-2} + \sigma^{-2}\sum_{i=1}^n \lambda_i^{-1} \ \ \ \mbox{and} \ \ \ \mu_1 = \tau_1^2\left(\tau_0^{-2}\mu_0+\sigma^{-2} \sum_{i=1}^n \lambda_i^{-1}y_i\right). \]

The draws \(\{(\mu,\lambda_1,\ldots,\lambda_n)^{(j)}\}_{j=1}^M\) are draws from \(p(\mu,\lambda_1,\ldots,\lambda_n|\mbox{data})\), while \(\{\mu^{(j)}\}_{j=1}^M\) are draws from \(p(\mu|\mbox{data})\).

4 Illustration

set.seed(31416)
mu.true = 0
sigma = 1
sigma2 = sigma^2
nu = 4
n = 30
y = mu.true+sigma*rt(n,df=nu)
ybar = mean(y)

mu0 = 0
tau02 = 2
mu = ybar

M = 10000
mu.d = rep(0,M)
for (iter in 1:M){
  lambda = 1/rgamma(n,(nu+1)/2,(nu+(y-mu)^2/sigma2)/2)
  tau12 = 1/(1/tau02+sum(1/lambda)/sigma2)
  mu1 = tau12*(mu0/tau02+sum(y/lambda)/sigma2)
  mu = rnorm(1,mu1,sqrt(tau12))
  mu.d[iter] = mu
}

par(mfrow=c(2,2))
plot(y,xlab="Observation",main="")
abline(h=ybar)
ts.plot(mu.d,xlab="iterations",main=expression(mu),ylab="")
acf(mu.d,main="")
hist(mu.d,prob=TRUE,main="",ylab="Posterior density",xlab=expression(mu))