Home > Bayesian computation > Normal approximation to posteriors

Normal approximation to posteriors

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 y_1, ..., y_{30} are independent where y_i is Poisson(\lambda_i), where the \lambda_i satisfy the log linear model

\log \lambda_i = \beta_0 + \beta_1 WEEKEND_i,

where WEEKEND_i is 0 (1) if the day is on a weekday (weekend).  I place a uniform prior on (\beta_0, \beta_1).

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)
contour
Advertisements
Categories: Bayesian computation
  1. No comments yet.
  1. No trackbacks yet.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: