digamma = function(x,a,b){dgamma(1/x,a,b)/(x^2)} prior = function(mu,sig){dnorm(mu,m0,sC0)*digamma(sig^2,a0,b0)} like = function(mu,sig){prod(dt((y-mu)/sig,df=nu)/sig)} y = c(-2.991,-2.845,-2.640,-1.962,-1.646,-1.627,-1.185,-0.997,-0.917,-0.726,-0.630, -0.382,-0.335,-0.300,-0.216,-0.178,-0.136,-0.013,0.185,0.255,0.367,0.422,0.440, 0.447,0.454,0.539,0.558,0.568,0.587,0.696,0.716,0.832,0.861,0.869,0.898,1.012, 1.073,1.081,1.112,1.159,1.332,1.342,1.402,1.560,1.728,1.785,2.364,2.491,2.558,3.575) n = length(y) # Gaussian data (simulated) with two extreme outliers # Data demeaned (sample mean zero) and scaled (sample std deviation one) set.seed(12345) n = 50 y = rnorm(n) y[y==min(y)]=-10 y[y==max(y)]=10 y = y-mean(y) y = y/sqrt(var(y)) ############################################################################ # Case 1: Fixing nu ############################################################################ nu = 5 # Prior hyperparameters m0 = 0 C0 = 1 sC0 = sqrt(C0) a0 = 5 b0 = 1 N = 200 mus = seq(-2.5,2.5,length=N) sigs = seq(0.001,1.5,length=N) likes = matrix(0,N,N) priors = matrix(0,N,N) posts = matrix(0,N,N) for (i in 1:N) for (j in 1:N){ likes[i,j] = like(mus[i],sigs[j]) priors[i,j] = prior(mus[i],sigs[j]) posts[i,j] = likes[i,j]*priors[i,j] } par(mfrow=c(1,1)) contour(mus,sigs,priors,xlab=expression(mu),ylab=expression(sigma),col=2, drawlabels=FALSE,lwd=2) contour(mus,sigs,likes,add=TRUE,drawlabels=FALSE,lwd=2,col=3) contour(mus,sigs,posts,add=TRUE,col=4,drawlabels=FALSE,lwd=2) legend("topleft",legend=c("Prior","Likelihood","Posterior"),col=2:4,lwd=2) # Gibbs sampler lambda = rep(1,n) sig2 = 1 M0 = 10000 M = 10000 niter = M0+M draws = matrix(0,niter,2) for (iter in 1:niter){ V = 1/(1/C0+sum(1/lambda)/sig2) m = V*(m0/C0+sum(y/lambda)/sig2) mu = rnorm(1,m,sqrt(V)) sig2 = 1/rgamma(1,a0+n/2,b0+sum((y-mu)^2/lambda)/2) lambda = 1/rgamma(n,(nu+1)/2,(nu+(y-mu)^2/sig2)/2) draws[iter,] = c(mu,sqrt(sig2)) } draws = draws[(M0+1):niter,] par(mfrow=c(2,3)) ts.plot(draws[,1],xlab="Iterations",ylab="",main=expression(mu)) acf(draws[,1],main="") hist(draws[,1],prob=TRUE,xlab="",main="") ts.plot(draws[,2],xlab="Iterations",ylab="",main=expression(sigma)) acf(draws[,2],main="") hist(draws[,2],prob=TRUE,xlab="",main="") par(mfrow=c(1,1)) plot(draws,xlab=expression(mu),ylab=expression(sigma)) contour(mus,sigs,posts,add=TRUE,col=6,drawlabels=FALSE,lwd=3) title("Joint posterior") draws0 = draws ############################################################################ # Case 2: Learning nu ############################################################################ # Gibbs sampler (including nu) nu.min = 1 nu.max = 100 c = 10 lambda = rep(1,n) sig2 = 1 nu = 5 M0 = 10000 M = 10000 niter = M0+M draws = matrix(0,niter,3) for (iter in 1:niter){ V = 1/(1/C0+sum(1/lambda)/sig2) m = V*(m0/C0+sum(y/lambda)/sig2) mu = rnorm(1,m,sqrt(V)) sig2 = 1/rgamma(1,a0+n/2,b0+sum((y-mu)^2/lambda)/2) lambda = 1/rgamma(n,(nu+1)/2,(nu+(y-mu)^2/sig2)/2) l1 = max(nu.min,nu-c) l2 = min(nu.max,nu+c) nu1 = sample(l1:l2,size=1) accept = min(0,dgamma(nu1,2,0.1,log=TRUE)+sum(dt((y-mu)/sqrt(sig2),df=nu1,log=TRUE))- dgamma(nu,2,0.1,log=TRUE)-sum(dt((y-mu)/sqrt(sig2),df=nu,log=TRUE))) if (log(runif(1))4 stdev = draws[ind,2]*sqrt(draws[ind,3]/(draws[ind,3]-2)) plot(draws[ind,3],stdev,xlab=expression(tau),ylab="Standard deviation",pch=16) par(mfrow=c(1,1)) xxx = seq(min(y),max(y),length=1000) plot(xxx,dnorm(xxx),type="l") lines(xxx,dt(xxx,df=median(draws[,3])),col=2)