# Mechanics for combining likelihood and prior in non-trivial case

Bayes rule is simple enough on its face:

$$pr(B|A) = \frac{pr(A|B)pr(B)}{pr(A)}$$

If these things are known scalar probabilities, the answer is simple to compute.

But I'm failing to understand what happens when I'm working with actual probability density functions.

For example, what if my prior ($$pr(B)$$) is $$\mathcal{N}\left(\left[\begin{array}{c}{ 1\\ 1 }\end{array}\right], \left[\begin{array}{cc} 1 & .5 \\ .5 & 1 \end{array}\right] \right)$$

Then I collect some new data that looks like this:

library(MASS)
set.seed(8675309)
N <- 1000

X <- mvrnorm(N, c(0, 0), matrix(c(1, .25, .25, 10), 2))
y <- X %*% c(.5, 2.5) + rnorm(N, sd = 25)

m <- lm(y~X-1)


The likelihood is also normal, with the mean being the coef vector

> coef(m)
X1        X2
0.4975087 2.5471237


and the variance being

> vcov(m) *sigma(m)
X1         X2
X1 16.4070280 -0.5447927
X2 -0.5447927  1.4990280


How is one to plug this into Bayes rule? This is what I've been struggling to understand.

Do I simply sample from each of these bivariate normals, multiplying and dividing? Do I have to worry about the fact that thir elements are correlated? In other words, am I just doing a hadamard product for multiplying the prior by the likelihood? And then the same thing in inverse to scale?

A simple simulation makes it seem like that's not the case; these results are weird:

AB <- mvrnorm(N, coef(m), vcov(m) *sigma(m))
B <- mvrnorm(N, c(1,1), matrix(c(1,.5,.5,1),2))
A <- mvrnorm(N, c(0, 0), matrix(c(1, .25, .25, 10), 2))

bhat <- AB*B/A

hist(bhat[,1])
hist(bhat[,2])


And where does MCMC come in? Is that just a way to figure out where in the distributions I should be sampling from?

Would appreciate an answer that could help me fill in the gaps!

EDIT I think that what I needed was for someone to tell me "do it over a grid" or "do it over the domain of your data." And "multiply the densities, point by point, not the values"

Here's a univariate example:

G <- seq(-5, 5, .01)
pB <- dnorm(x = G, mean = -3, sd = .5)
pB <- pB/sum(pB)
plot(G, pB, type = "l")
pAB <- dnorm(x = G, mean = 3, sd = .5)
pAB <- pAB/sum(pAB)
lines(G, pAB, type = "l", col = "red")
pBA <- pB*pAB/(sum(pB*pAB))
lines(G, pBA, type = "l", col = "blue")


But I must be doing something wrong here: why is the posterior pointier than the prior and the likelihood in this case? Or is it that the normalization something more than a factor to scale the PDF to sum to 1?

But I must be doing something wrong here: why is the posterior pointier than the prior and the likelihood in this case? Or is it that the normalization something more than a factor to scale the PDF to sum to 1?

• Remind me in an hour to answer this. Aug 5 '19 at 16:37
• @DemetriPananos, please and thank you. Aug 5 '19 at 18:18
• Regarding the posterior being "pointier" than the prior: this is what you want to see. The spread of the normal distribution for B reflects your uncertainty in its value, so, if the data gives any information, the posterior has a smaller spread than the prior, due to the information decreasing the uncertainty. If the spread doesn't change, the data gives no information. If the spread increases, something's gone wrong with the calculation, because information can't increase the uncertainty. Aug 5 '19 at 20:08
• @AccidentalStatistician but what if my likelihood has a higher variance than my prior? Intuitively, shouldn't the distribution widen with evidence against the notion of certainty? Aug 6 '19 at 1:13
• @user8675309 The posterior would be wider for a wider likelihood than for a narrower likelihood, but it shouldn't be wider than the prior. Take the extreme case, where the likelihood has infinite variance. This likelihood is flat, i.e. any parameter value has equal likelihood. The data therefore gives no information about the parameter. The posterior is therefore the same as the prior. Since this is the extreme case, the posterior therefore can't be wider than the prior. Aug 6 '19 at 16:18

I think your problem setup is a little troublesome. Let's take a step back and write down a model.

Bayes' Rule says

$$p(\theta \vert y) \propto p(y \vert \theta) p(\theta)$$

Equivalently

$$\log(p(\theta \vert y)) \propto \log(p(y \vert \theta)) + \log( p(\theta))$$

I'm going to compute things on the log scale because it's easier.

Let's say I have some data that I think comes from a bivariate normal with known covariance matrix (assume it is the identity for simplicity). Our job is to estimate the mean of the bivariate normal from data. Our likelihood is then

$$y \vert \mu \sim \mathcal{N}(\mu, \Sigma)$$

Now we need a prior for mu. Let's just use the one you gave

$$\mu \sim \mathcal{N}\left(\left[\begin{array}{l}{1} \\ {1}\end{array}\right],\left[\begin{array}{cc}{1} & {.5} \\ {.5} & {1}\end{array}\right]\right)$$

Now, we need to evaluate the posterior along a 2d grid of possible values of $$\mu$$. To compute the posterior probability (well, not really, just something that is proportional to the probability) we compute

$$\log(p(y \vert \mu)) + \log( p(\mu))$$

But the important part is that we need to compute this over a grid. That grid will be the elements of $$\mu$$. Let's turn to python to do this

import numpy as np
from scipy.stats import multivariate_normal
import matplotlib.pyplot as plt
#Data Generating Parameters
true_mean = [-1.5, -1.75]
true_Sigma = np.eye(2)

#Generate data from our data generating process
data = multivariate_normal.rvs(mean = true_mean, cov = true_Sigma, size = 100)

#Our prior
prior = multivariate_normal(mean = np.ones(2), cov = np.array([[1,0.5], [0.5,1]]))

#Generate a grid on which to evaluate the posterior
grid = np.linspace(-2,2,26)
mu_1, mu_2 = np.meshgrid(grid,grid)
#Reshape the grid for ease of computation
MU = np.stack((mu_1, mu_2), axis = -1).reshape(-1,1,2)

#Instantiate the loglikelihood
loglik = prior.logpdf(MU)

fig,ax = plt.subplots(dpi = 120)
ax.axes.set_aspect('equal')
plt.contourf(mu_1, mu_2,loglik.reshape(mu_1.shape))
plt.title('Prior')
plt.show()

LL = []
for m in MU:
#Evualuate the likelihood of the data over the grid
LL.append(multivariate_normal.logpdf(data, mean = m[0], cov = np.eye(2)))

LL = np.array(LL)
#Add the log likelihood to the log prior
loglik+=LL.sum(axis = 1)

fig,ax = plt.subplots(dpi = 120)
ax.axes.set_aspect('equal')
plt.contourf(mu_1, mu_2,loglik.reshape(mu_1.shape))
plt.title('Posterior')
plt.show()


This code produces the following

As you can see, the prior looks like your prior distribution (as it should). After seeing data, the posterior has shifted the density towards the true mean of the data generating process. So it isn't enough to just evaluate bayes rule as you do in your samples. You have to work with densities, which can get tricky.

And where does MCMC come in? Is that just a way to figure out where in the distributions I should be sampling from?

As you can imagine, the densities become very difficult to deal with with a modest amount of parameters. MCMC and HMC (Hamiltonian Monte Carlo) are ways to generating markov chains with limiting distributions equal to our posterior. There is a lot to say about MCMC and HMC, so I suggest you create another question if you are really interested.

• Thanks, that's helpful. What if I want a true PDF, rather than something that is just proportional to the posterior? Do I just divide by the sum of the densities, after exponentiating them from the log scale? Aug 5 '19 at 19:44
• @user8675309 You won't be able to obtain the entire density since the normal is supported on the reals. You can either a) extend the grid to a large number and then normalize the loglik to sum to 1 after exponentiating, thereby obtaining an approximation to the posterior, or b) do compute the analytical posterior (which shouldn't be that hard. See chapter 3 of Bayesian Data Analysis). Aug 5 '19 at 20:18