# MC Sampling from a Posterior Distribution

Bayes rule is:

$p(x|D)\propto p(D|x)p(x)$

Imagine $p(x)$ and $p(D|x)$ are both Gaussian distributed. I understand in this case that I can obtain the posterior analytically (conjugate priors). But I'm a bit confused on how to sample from the posterior distribution via an MC (Monte Carlo) method.

If $p(x)$ is Gaussian, I can easily pull a sample from there, but then do I put this sample into the likelihood function? If that's the case I'm evaluating a probability value at a particular point, and not actually "taking a sample". I'm a little confused on this front. Basically once I take a sample from $p(x)$ what is the procedure I need to do with that sample (i.e. where / how do I pass it on).

I would just like to see a histogram of this posterior distribution.

• To the best of my knowledge, you can either choose a proposal distribution and do Metropolis-Hastings or sample from the prior and follow the method of the approximate Bayesian computation(ABC). But if there is only one parameter that needs to be sample, why not just use rejection sampling or importance sampling method? – Daeyoung Lim Sep 6 '16 at 10:31
• Rejection sampling sounds like a good approach. However I was thinking of hopefully a way in which things work more like a graphical model. So for example if we have $p(x|a)p(a|b)p(b|c)p(c)$ can we sample from the $p(c)$ prior and then pass that through all the other probability measurements sequentially and obtain one final posterior sample? Rejection sampling would allow for sampling, but I feel there should be some way to do it in this graphical model method. :/ – pche8701 Sep 6 '16 at 14:04
• Actually I completely missed on that approximate Bayesian computation (ABC) method you mentioned. That sounds like what I may be looking for? In your experience will this allow for sampling in the way I described in my previous comment? – pche8701 Sep 6 '16 at 14:06
• It seems like there are more than one parameter. Then you should think of MCMC algorithms like Gibbs samplers or Metropolis-Hastings updates. – Daeyoung Lim Sep 6 '16 at 14:07
• actually, if you have a lot of parameters, as is the case with graphical models, you're subject to the curse of dimensionality. So you wouldn't want to do ABC in such cases. – Daeyoung Lim Sep 6 '16 at 14:08

You worte

If $$p(X)$$ is Gaussian, I can easily pull a sample from there, but then do I put this sample into the likelihood function? If that's the case I'm evaluating a probability value at a particular point, and not actually "taking a sample

I think you misunderstand the concept of conditional probability, when we sample from the prior and plug the sampled value into the likelihood , we do not evaluate a probability value at a particular point we just get a conditional distribution given the sampled value from the prior. Suppose that the target distribution is the joint probability distribution of two variables, a discrete variable $$\delta \in \{1,2,3\}$$ and a continuous variable $$\theta \in \mathbb{R}$$ then the target density will be defined as :

01-The discrete part of this density is : $$\{Pr(\delta =1), Pr(\delta = 2), Pr( \delta= 3)\} = (.45, .10, .45)$$

02- The continuous part of this density is $$P(\theta|\delta) = dnorm(\theta,\mu_{\delta},\sigma_{\delta })$$, where $$(\mu_1,\mu_2,\mu_3) = (−3, 0, 3)$$ and $$(\sigma^2_1,\sigma^2_2,\sigma^2_3) = (1/3, 1/3, 1/3)$$.

This is a mixture of three normal densities and a plot of the exact marginal density of ,$$p(\theta) =\sum p(\theta|\delta)p(\delta)$$ will be like that

that can be done using the following code in R :

mu<-c(-3,0,3)
s2<-c(.33,.33,.33)
w<-c(.45,.1,.45)

ths<-seq(-5,5,length=100)
plot(ths, w[1]*dnorm(ths,mu[1],sqrt(s2[1])) +
w[2]*dnorm(ths,mu[2],sqrt(s2[2])) +
w[3]*dnorm(ths,mu[3],sqrt(s2[3])) ,type="l" )


Now it is very easy to obtain independent Monte Carlo samples from the joint distribution of $$(\delta,\theta)$$ First, a value of $$\delta$$ is sampled from its marginal distribution, then the value is plugged into $$p(\theta|\delta)$$ from which a value of $$\theta$$ is sampled (Note the sampled pair $$(\delta,\theta)$$ represents a sample from the joint distribution $$p(\delta,\theta)=p(\delta)p(\theta|\delta)$$) then the empirical distribution of the $$\theta$$ samples provides an approximation to the marginal distribution $$p(\theta) =\sum p(\theta|\delta)p(\delta)$$ . Now a histogram of 1,000 Monte Carlo -values generated in this way or in other words the empirical distribution of the Monte Carlo samples looks like $$p(\theta)$$.

set.seed(1)
S<-2000
d<-sample(1:3,S, prob=w,replace=TRUE)
th<-rnorm(S,mu[d],sqrt(s2[d]))
THD.MC<-cbind(th,d)

ths<-seq(-6,6,length=1000)
plot(ths, w[1]*dnorm(ths,mu[1],sqrt(s2[1])) +
w[2]*dnorm(ths,mu[2],sqrt(s2[2])) +
w[3]*dnorm(ths,mu[3],sqrt(s2[3])) ,type="l" , xlab=expression(theta),ylab=
expression( paste( italic("p("),theta,")",sep="") ),lwd=2 ,ylim=c(0,.40))
lines( ths, w[1]*dnorm(ths,mu[1],sqrt(s2[1])) +
w[2]*dnorm(ths,mu[2],sqrt(s2[2])) +
w[3]*dnorm(ths,mu[3],sqrt(s2[3])),lwd=2 )


For more details see A First Course in Bayesian Statistical Methods the book provides also a very good comparison between MC and MCMC Chapter 6 section "Introduction to MCMC diagnostics ".

• That is a very good way to look at it! The joint between a discrete and a continuous variable. I'm not sure how this is frequentist though, it seems very Bayesian (prior is $p(\delta)$, likelihood is $p(\theta|\delta)$)? – pche8701 Sep 8 '16 at 23:59
• @pche8701 see edit – Bahgat Nassour Sep 10 '16 at 7:25

Note that you have a sample, i.e. the data, and you want to find the credibility of an unknown parameter, below I explain a special kind of MC nemely MCMC, markov Chain Monte Carlo.

Maybe it's good to make the notation a bit more explicit, and for simplicity take the one-dimensional case (this is only for simplicity, it also works in the multi-dimensional case). You say that your likelihood is Gaussian just as well as your prior. If we take the one-dimensional case then our goal it to estimate one parameter, e.g. the unkown mean $\mu$.

The posterior is then $P(\mu|_D) \propto P(D|_\mu) p(\mu)$

Now in MCMC (Monte Carlo Markov Chain), you start from an initial value $\mu=x_1$ and make a proposal move $x_2$, in order to decide whether the proposal will be accepted, you should evaluate the posterior at $x_1$ and at $x_2$.

So you have to compute $P(D|_\mu=x_1) p(\mu=x_1)$ (and the same for $x_2$). The prior is known (you have chosen it yourself), so you can evaluate it at $x_1$.

The likelihood $P(D|_\mu=x_1)$ at $x_1$ is also knwown; you know it is gaussian and (in the one dimensional case) only the mean is unknown but you say that you want to evaluate it at $\mu=x_1$ so you have the mean and you can compute the value for $P(D|_\mu=x_1)$ (the data $D$ is also known of course). Similar for $X_2$. So it is the 'Gaussian formula'' where you plug in $\mu=x_1$ and $X=D$ ($\sigma$ is known because we assumed the one dimensional case for simplicity, so only one parameter, $\mu$ is unknown, but you evaluate at $\mu=x_1$, so that's solved).

So we can compute (up to a constant) $\pi_1=P(\mu=x_1|_D)$ and $\pi_2=P(\mu=x_2|_D)$.

Remember that $x_2$ was only a proposal move. This proposal is accepted if $\pi_2$ > $\pi_1$ , if $\pi_2 \le \pi_1$ then the move is accepted with a probablity $\frac{\pi_2}{\pi_1}$ (P.S. a symmetric transition kernel is assumed).

If the move is accepted then move to $x_2$ else stay at $x_1$ and next generate the next proposal, ...

This algorithm is constructed such that the so-called ''balance equations'' are fulfilled, and that is suffcient for the chain to converge to the posterior ''in the end''.

Summary: if you can't sample from the posterior directly, then construct a markov chain that has the posterior as limiting distribution and use this Markov chain to ' walk to' that limiting distribution. The algorithm supra defines such a markov chain using a special case of the Metropolis-Hastings procedure.

[[Note, if your data is a sample with multiple elements, then you take a product of likelihoods, assuming independence). ]]