Linear regression with Student’s t errors

We collected data from 15 volunteers (my students and myself) and below we summarize the data information graphically and entertain two competing simple linear regression models: i) Gaussian errors, ii) Student’s \(t\) errors.

hedi = 115
n    = 15
data = matrix(c(178,70,175,80,184,82,200,89,183,80,173,65,185,70,
162,70,190,85,174,80,180,87,165,74,176,72,174,63,196,hedi),n,2,byrow=TRUE)

x = data[,1]
y = data[,2]

loglike.t = function(alpha,beta,sig,nu){
    tau = sqrt((nu-2)*sig^2/nu)
  sum(dt((y-alpha-beta*x)/tau,df=nu,log=TRUE))-n*log(tau)
}

prior.nu = function(nu){
  sqrt((nu+1)/(nu+3)*nu/(nu+3)*(
    trigamma(nu/2)-trigamma((nu+1)/2)-2*(nu+3)/(nu*(nu+1)^2)))
}

nus = 3:30
priornu = prior.nu(nus)
priornu = priornu/sum(priornu)

M       = 1000000
alpha.d = runif(M,-250,150)
beta.d  = runif(M,-0.5,2)
sig.d   = runif(M,1,30)
nu      = sample(nus,size=M,replace=TRUE,prob=priornu)
w       = rep(0,M)
for (i in 1:M)
  w[i] = loglike.t(alpha.d[i],beta.d[i],sig.d[i],nu[i])
w   = exp(w-max(w))
w   = w/sum(w)
ind = sample(1:M,size=M,replace=TRUE,prob=w)
pars = cbind(alpha.d[ind],beta.d[ind],sig.d[ind],nu[ind])


tab = t(apply(pars,2,quantile,c(0.05,0.25,0.5,0.75,0.95)))
row.names(tab) = c("alpha","beta","sigma","nu")
tab

##                 5%         25%         50%         75%       95%
## alpha -129.8864379 -73.8993458 -39.2038016 -11.9204183 29.222764
## beta     0.2703433   0.5018516   0.6514483   0.8431003  1.158998
## sigma    7.4547320   9.3565316  11.0873509  13.4132851 18.255269
## nu       3.0000000   3.0000000   5.0000000   9.0000000 21.000000

par(mfrow=c(2,2))
hist(pars[,1],prob=TRUE,main=expression(alpha),xlab="")
hist(pars[,2],prob=TRUE,main=expression(beta),xlab="")
hist(pars[,3],prob=TRUE,main=expression(sigma),xlab="")
#hist(pars[,4],prob=TRUE,main=expression(nu))
plot(table(pars[,4])/M,main=expression(nu),ylab="Probability")
lines(nus+0.25,priornu,type="h",col=2,lwd=2)
legend("topright",legend=c("Prior","Posterior"),col=2:1,lty=1,lwd=2,bty="n")