############################################################## # Sampling Importance Resampling (SIR) ############################################################## # # We model the mean of Gaussian data via # a Laplace (Double Exponential) prior, ie. # # y1,...,yn ~ N(mu,1) # # mu ~ DE(b,B) # ############################################################## rm(list=ls()) DE = function(x,b,B){ exp(-abs(x-b)/B)/(2*B) } # Simulated data set.seed(15497) n = 20 beta = 0 y = rnorm(n,beta,1) # Laplace prior b = 0 B = 2 # Proposal 1: normal prior c0 = 0 D0 = 900 M = 1000000 betas0 = rnorm(M,c0,sqrt(D0)) # Proposal 2: normal posterior based on a normal prior # but inflating its variance by a factor of 100 D1 = 1/(1/D0+n) c1 = D1*(c0/D0+sum(y)) D1 = 100*D1 betas = rnorm(M,c1,sqrt(D1)) # Computing weights w0 = rep(0,M) w = rep(0,M) for (i in 1:M){ w0[i] = DE(betas0[i],b,B)*exp(-0.5*sum((y-betas0[i])^2))/dnorm(betas0[i],c0,sqrt(D0)) w[i] = DE(betas[i],b,B)*exp(-0.5*sum((y-betas[i])^2))/dnorm(betas[i],c1,sqrt(D1)) } w0 = w0/sum(w0) w = w/sum(w) # Resampling N = 0.1*M ind0 = sample(M,size=N,prob=w0,replace=TRUE) ind = sample(M,size=N,prob=w,replace=TRUE) betas0.draw = betas0[ind0] betas.draw = betas[ind] L = min(betas0.draw,betas.draw) U = max(betas0.draw,betas.draw) pdf(file="workedexample-posterior-via-SIR-using-Laplaceprior.pdf",width=8,height=8) par(mfrow=c(2,2)) hist(betas0,prob=TRUE,main=paste("Proposal 1: N(",c0,",",D0,")",sep=""),xlab=expression(beta));box() hist(betas0.draw,xlim=c(L,U),prob=TRUE,main="Posterior - Proposal 1",xlab=expression(beta));box() hist(betas,prob=TRUE,main=paste("Proposal 2: N(",round(c1,2),",", round(D1,2),")",sep=""),xlab=expression(beta));box() hist(betas.draw,xlim=c(L,U),prob=TRUE,main="Posterior - Proposal 2",xlab=expression(beta));box() dev.off()