How can I compute a posterior density estimate from a prior and likelihood?

Question

I am trying to understand how to use Bayes' theorem to calculate a posterior but am getting stuck with the computational approach, e.g., in the following case it is not clear to me how to take the product of the prior and likelihood and then calculate the posterior:

For this example, I am interested in calculating the posterior probability of $\mu$ and I use a standard normal prior on $\mu$ $p(\mu)\sim N(\mu = 0, \sigma = 1)$, but I want to know how to calculate the posterior from a prior on $\mu$ that is represented by an MCMC chain, so I will use 1000 samples as my starting point.

• sample 1000 from the prior.

  set.seed(0)
  prior.mu      <- 0
  prior.sigma   <- 1
  prior.samples <- sort(rnorm(1000, prior.mu, prior.sigma))
  

• make some observations:

  observations <- c(0.4, 0.5, 0.8, 0.1)
  

• and calculate the likelihood, e.g. $p(y | \mu, \sigma)$:

  likelihood <- prod(dnorm(observations, mean(prior.samplse), sd(prior.samples)))
  

what I don't quite understand is:

1. when / how to multiply the prior by the likelihood?
2. when / how to normalize the posterior density?

please note: I am interested in the general computational solution that could be generalizable problems with no analytical solution

Answer 1

You have several things mixed up. The theory talks about multiplying the prior distribution and the likelihood, not samples from the prior distribution. Also it is not clear what you have the prior of, is this a prior on the mean of something? or something else?

Then you have things reversed in the likelihood, your observations should be x with either prior draws or known fixed constants as the mean and standard deviation. And even then it would really be the product of 4 calls to dnorm with each of your observations as x and the same mean and standard deviation.

What is really no clear is what you are trying to do. What is your question? which parameters are you interested in? what prior(s) do you have on those parameters? are there other parameters? do you have priors or fixed values for those?

Trying to go about things the way you currently are will only confuse you more until you work out exactly what your question is and work from there.

Below is added due after the editing of the original question.

You are still missing some pieces, and probably not understanding everything, but we can start from where you are at.

I think you are confusing a few concepts. There is the likelihood that shows the relationship between the data and the parameters, you are using the normal which has 2 parameters, the mean and the standard deviation (or variance, or precision). Then there is the prior distributions on the parameters, you have specified a normal prior with mean 0 and sd 1, but that mean and standard deviation are completely different from the mean and standard deviation of the likelihood. To be complete you need to either know the likelihood SD or place a prior on the likelihood SD, for simplicity (but less real) I will assume we know the likelihood SD is $\frac12$ (no good reason other than it works and is different from 1).

So we can start similar to what you did and generate from the prior:

> obs <- c(0.4, 0.5, 0.8, 0.1)
> pri <- rnorm(10000, 0, 1)

Now we need to compute the likelihoods, this is based on the prior draws of the mean, the likelihood with the data, and the known value of the SD. The dnorm function will give us the likelihood of a single point, but we need to multiply together the values for each of the observations, here is a function to do that:

> likfun <- function(theta) {
+ sapply( theta, function(t) prod( dnorm(obs, t, 0.5) ) )
+ }

Now we can compute the likelihood for each draw from the prior for the mean

> tmp <- likfun(pri)

Now to get the posterior we need to do a new type of draw, one approach that is similar to rejection sampling is to sample from the prior mean draws proportional to the likelihood for each prior draw (this is the closest to the multiplication step you were asking about):

> post <- sample( pri, 100000, replace=TRUE, prob=tmp )

Now we can look at the results of the posterior draws:

> mean(post)
[1] 0.4205842
> sd(post)
[1] 0.2421079
> 
> hist(post)
> abline(v=mean(post), col='green')

and compare the above results to the closed form values from the theory

> (1/1^2*mean(pri) + length(obs)/0.5^2 * mean(obs))/( 1/1^2 + length(obs)/0.5^2 )
[1] 0.4233263
> sqrt(1/(1+4*4))
[1] 0.2425356

Not a bad approximation, but it would probably work better to use a built-in McMC tool to draw from the posterior. Most of these tools sample one point at a time not in batches like above.

More realistically we would not know the SD of the likelihood and would need a prior for that as well (often the prior on the variance is a $\chi^2$ or gamma), but then it is more complicated to compute (McMC comes in handy) and there is no closed form to compare with.

The general solution is to use existing tools for doing the McMC calculations such as WinBugs or OpenBugs (BRugs in R gives an interface between R and Bugs) or packages such LearnBayes in R.