rm(list=ls())

ldig = function(x,a,b){
  dgamma(1/x,a,b,log=TRUE)-2*log(x)
}

dt.hedi = function(x,mu,sig,nu){
  dt((x-mu)/sig,df=nu)/sig
}

ldt.hedi = function(x,mu,sig,nu){
  dt((x-mu)/sig,df=nu,log=TRUE)-log(sig)
}

lpost = function(y,x,beta,sig2,nu,b0,B0,nu0,sig20){
  ldig(sig2,nu0/2,nu0*sig20/2)+
  dnorm(beta,b0,sqrt(B0),log=TRUE)+
  sum(ldt.hedi(y,x*beta,sqrt(sig2),nu))
}

set.seed(12345)
n    = 100
beta = 2.0
sig2 = 1.0
nu   = 2
sig  = sqrt(sig2)
erro = sig*rt(n,df=nu)
x    = runif(n)
y    = beta*x + erro

par(mfrow=c(1,1))
plot(x,y)
abline(lm(y~x-1),col=2,lwd=2)


# hyperparameters of the priors
b0    = 0.0
B0    = 100
nu0   = 5
sig20 = 1.0
numax = 50

# initial values (for my MCMC scheme)
b  = 0
s2 = 1
d  = 10

# tuning standard deviation
stdev = 0.1

# Running the MCMC
set.seed(2468)
burn  = 10000
M     = 2000
LAG   = 50
niter = burn+LAG*M
betas = rep(0,niter)
sig2s = rep(0,niter)
nus   = rep(0,niter)
for (iter in 1:niter){
  b1 = rnorm(1,b,stdev)
  A  = lpost(y,x,b1,s2,nu,b0,B0,nu0,sig20)-lpost(y,x,b,s2,nu,b0,B0,nu0,sig20)
  logacceptance = min(0,A)
  if (log(runif(1))<logacceptance){
    b = b1
  }
  betas[iter] = b
  s21 = rgamma(1,1,1)
  A  = lpost(y,x,b,s21,nu,b0,B0,nu0,sig20)-lpost(y,x,b,s2,nu,b0,B0,nu0,sig20)
  B  = dgamma(s2,1,1,log=TRUE)-dgamma(s21,1,1,log=TRUE)
  logacceptance = min(0,A+B)
  if (log(runif(1))<logacceptance){
    s2 = s21
  }
  sig2s[iter] = s2
  d1 = sample(1:numax,size=1,prob=rep(1/numax,numax))
  A  = lpost(y,x,b,s21,d1,b0,B0,nu0,sig20)-lpost(y,x,b,s2,d,b0,B0,nu0,sig20)
  logacceptance = min(0,A)
  if (log(runif(1))<logacceptance){
    d = d1
  }
  nus[iter] = d
}
betas = betas[seq(burn+1,niter,by=LAG)]
sig2s = sig2s[seq(burn+1,niter,by=LAG)]
nus   = nus[seq(burn+1,niter,by=LAG)]

pdf(file="regression-with-t-error-inclass.pdf",width=10,height=10)
par(mfrow=c(3,3))
ts.plot(betas,ylab="",main=expression(beta))
acf(betas,main="")
hist(betas,prob=TRUE,xlab="",main="");box()
ts.plot(sig2s,ylab="",main=expression(sigma^2))
acf(sig2s,main="")
hist(sig2s,prob=TRUE,xlab="",main="");box()
ts.plot(nus,ylab="",main=expression(nu))
acf(nus,main="")
plot(table(nus)/M,ylab="Probability",xlab="")
dev.off()


mbeta = mean(betas)
xx    = seq(min(x),max(x),length=100)
fit   = xx*mbeta
par(mfrow=c(1,1))
plot(x,y)
abline(lm(y~x-1),col=2,lwd=2)
lines(xx,fit,col=4,lwd=2)