Student’s t model

Student’s \(t\) model

Assume that \(y_1,\ldots,y_n\) are iid \(t_\nu(0,\sigma^2)\), for known number of degrees of freedom \(\nu\), and that the prior of \(\sigma^2\) is inverse-gamma with parameters \(a\) and \(b\), i.e. \(\sigma^2 \sim IG(a,b)\). Therefore, \[ p(\sigma^2|y_1,\ldots,y_n,\nu) \propto p(\sigma^2|a,b)\prod_{t=1}^n p(y_i|\sigma^2,\nu), \] is of unknown form.

Simulating some artificial data

set.seed(2325)
n   = 100
sig = 2
nu  = 4
y   = sig*rt(n,df=nu)
boxplot(y,horizontal=TRUE)

SIR to draw from \(p(\sigma^2|y_1,\ldots,y_n,\nu)\)

Below we implement SIR by sampling from a highly inefficient candidate density, \(U(0,100)\).

dt.hedi = function(sig2){
  prod(dt(y/sqrt(sig2),df=nu)/sqrt(sig2))
}
set.seed(4321)
M = 10000
sig2.t = runif(M,0,100)
w = rep(0,M)
for (i in 1:M)
  w[i] = dt.hedi(sig2.t[i])
sig2 = sample(sig2.t,size=M,replace=T,prob=w)
hist(sig2,xlab=expression(sigma^2),prob=TRUE)
abline(v=sig^2,col=2,lwd=3)

mean(sig2)

## [1] 4.738057

Computing marginal likelihoods (aka normalizing constants, aka prior predictives)

Let us know compute \[ p(y_1,...,y_n|\nu) = \int p(y_1,\ldots,y_n|\sigma^2,\nu)p(\sigma^2|a,b)d\sigma^2 \] for \(\nu \in \{1,\ldots,k\}\) for a large \(k\), say \(k=100\). We will approximate the integral by simple Monte Carlo Integration by noticing that \[ p(y_1,\ldots,y_n|\nu) \approx \frac{1}{M} \sum_{i=1}^M p(y_1,\ldots,y_n|\sigma^{2(i)},\nu), \] where \(\{\sigma^{2(i)}\}_{i=1}^M\) are draws from \(p(\sigma^2|a,b)\), the prior for \(\sigma^2\) with \(a=b=3/2\).

dt.hedi = function(sig2,nu){
  sum(dt(y/sqrt(sig2),df=nu,log=TRUE)-0.5*log(sig2))
}

M = 1000
sig2 = 1/rgamma(M,3.5,3.5)
nu.max = 100
nus = 1:nu.max
logpred = matrix(0,M,nu.max)
for (j in 1:M)
for (i in 1:nu.max)
  logpred[j,i] = dt.hedi(sig2[j],nus[i])
A = max(logpred)
logpred1 = logpred-A
pred1 = exp(logpred1)
logpred1 = log(apply(pred1,2,sum))+ A -log(M)
plot(nus,logpred1,ylab="Log predictive",xlab=expression(nu))
abline(v=nu,col=2,lwd=3)