Bayesian Thinking

Normal approximation to posteriors

Advertisements

Today we talked about normal approximations to posterior distributions.  As an example, I observe the daily hits on my book’s website during each of the days one recent month.  I collect the following counts for the weekdays (Monday through Friday)

20,30,22,20,20,17,21,26,22,30,36,15,30,27,22,23,18,
24,28,23,12

20,30,22,20,20,17,21,26,22,30,36,15,30,27,22,23,18,24,28,23,12

and I collect the following counts for the weekends (Saturday and Sunday)

7,12,11,12,12,17,18,20,17

We assume that the counts are independent where is Poisson(), where the satisfy the log linear model

,

where is 0 (1) if the day is on a weekday (weekend).  I place a uniform prior on .

I program the log posterior in the function loglinearpost — the data matrix has two columns, where the first column is the values of WEEKEND and the second column has the counts.

loglinearpost=function (beta, data)
{
x = data[, 1]
y = data[, 2]
logf = function(x,  y, beta) {
lp = beta[1] + beta[2] * x
dpois(y,exp(lp),log=TRUE)}
sum(logf(x,y, beta))
}

I put the counts and the weekend indicators in the vectors y and x and form the data matrix.
y=c(7,12,11,12,12,17,18,20,1720,30,22,20,20,17,21,26,22,30,36,15,30,27,22,23,18,
24,28,23,12)
x=c(rep(0,9),rep(1,21))
data=cbind(x,y)
I obtain the normal approximation using the laplace function.  The output is the posterior mode and the corresponding variance-covariance matrix.
fit=laplace(loglinearpost,c(0,0),data)
fit
$mode
[1] 2.6390946 0.5025815
$var
[,1]         [,2]
[1,]  0.00793621 -0.007936210
[2,] -0.00793621  0.009993843
Below using the mycontour function, I show the exact posterior in RED and the approximating normal density in BLUE.  Note that I define a new function mydnorm that computes the log of the multivariate normal density (defined in dmnorm).   Note that the two sets of contours are very close, indicating that the normal approximation is very accurate in this case.
mycontour(loglinearpost,c(2.2,3,.1,.9),data,col=”red”,
xlab=”BETA0″,ylab=”BETA1″)
mydnorm=function(y,d)
dmnorm(y,d$mu,d$sigma,log=TRUE)
mycontour(mydnorm,c(2.2,3,.1,.9),
list(mu=fit$mode,sigma=fit$var),col=”blue”,add=TRUE)
Advertisements