Product of two normal distributions (for Bayes Rule) is not product of normal output variables?

Question

When we apply Bayes' rule in machine learning, we want to compute the posterior probability $P(y|X)$ by multiplying two probability distributions (the observed class-conditional likelihood $P(X|y)$ and the prior probability $P(y)$) and then dividing by the evidence $P(X)$, like so:

$$ P(y|X) = \frac{P(X|y) \times P(y)}{P(X)} $$

Suppose that the likelihood and prior terms are both Gaussian normal. Let's say $P(X|y) = \mathcal{N}(10.0, 4.0)$ and $P(y) = \mathcal{N}(12.0, 4.0)$. So we need to compute the product of these two Gaussian distributions.

Ater watching an exciting video

I've learned that the product of two Gaussian distributions is another Gaussian::

$$ \mathcal{N}(\mu_{1}, \sigma_1^{2}) \times \mathcal{N}(\mu_{2}, \sigma_1^{2}) = \mathcal{N}(\mu_{3}, \sigma_3^{2}) $$ where

\begin{align} \mu_3 &= \frac{\sigma_2^{2} \mu_1 + \sigma_1^{1} \mu_2}{\sigma_1^{2} + \sigma_2^{2}} \\ \sigma_3^{2} &= \frac{1}{\frac{1}{\sigma_1^{2}} + \frac{1}{\sigma_2^{2}} } \end{align}

In this example, $\mathcal{N}(10.0, 4.0) \times \mathcal{N}(12.0, 4.0) = \mathcal{N}(11.0, 2.0)$

However, that last formula is not what I get when I multiply the two Gaussian series 1 and 2. Here's some Python code to demonstrate what I mean.

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import math

def normal(x, mu, sigma_sq):
    # Computes Gaussian probability from
    # https://en.wikipedia.org/wiki/Gaussian_function

    normalizer = 1.0 / np.sqrt(2.0 * math.pi * sigma_sq)
    exp_term = - 0.5 * (x - mu)**2 / sigma_sq
    return normalizer * np.exp(exp_term)

x = np.arange(0, 20, 0.2)
y1 = normal(x, 10, 4)
y2 = normal(x, 12, 4)
y3 = normal(x, 11, 2)
y4 = y1 * y2

So:

• y1 is $\mathcal{N}(\mu_{1}, \sigma_1^{2})$

• y2 is $\mathcal{N}(\mu_{2}, \sigma_1^{2})$.

• y3 is $\mathcal{N}(\mu_{1}, \sigma_1^{2}) \times \mathcal{N}(\mu_{2}, \sigma_1^{2}) = \mathcal{N}(\mu_{3}, \sigma_3^{2})$.

• And y4 is the element-wise product y1 * y2.

Now, when I plot these four series y1 (blue), y2 (yellow), y3 (green), and y4 (red), the results of y3 and y4 are not the same!

plt.figure(figsize=(10, 7))
plt.plot(x, y1, label='y1 = N(10, 4)', color='b')
plt.plot(x, y2, label='y2 = N(12, 4)', color='y')
plt.plot(x, y3, label='y3 = N(10, 4) x N(12, 4) = N(11, 2)', color='g')
plt.plot(x, y4, label='y4 = y1 * y2', color='r')
plt.legend()
plt.grid(True)

Shouldn't y3 and y4 be the same? Why aren't they?

Answer 1

The problem is that you’re abusing the notation. What do you exactly mean by $\mathcal{N}(\mu, \sigma^2)$ in here? Is it values of random variable, it’s probability density function, or maybe cumulative distribution function?

You can multiply the probability density functions, this would be equal to their joint density if they are independent

$$ p_X(x) p_Y(y) = p_{XY}(x, y) $$

this is probability of observing $x$ and $y$ happening together. If there would be dependence between them, the results won’t hold by definition. It’s not the same as observing some value $z=xy$, notice that $z$ can have many divisors, so there could be multiple values of $x$ and $y$ that can lead to observing it, so obviously you won’t calculate the probability from probabilities of two particular divisors alone.

For example, imagine that you have two discrete random variables that both can take values in $\{1, 2, 3, 4\}$, $P(XY=4)$ is a probability that depends on observing any of the three pairs of $x$ and $y$ values $(1,4)$, $(2,2)$, $(4,1)$, so it would be equal to

$$P(X=1,Y=4) + P(X=2,Y=2) + P(X=4,Y=1)$$

Notice that to calculate $P(X\times Y)$ we summed the probabilities, math on the values does not apply to the math on the distributions.

Your Bayes theorem example is about multiplying conditional probability density $p(x|y)$ by marginal probability density $p(y)$ and dividing by marginal probability density $p(x)$ and the result is another conditional probability density $p(y|x)$. It’s not just multiplying two Gaussian pdf’s to get a third one. Your example didn’t work because you didn’t apply Bayes theorem correctly. Obviously

$$p(x|y)p(y) \ne \frac{p(x|y)p(y)}{p(x)}$$

Moreover, the fact that the result is a Gaussian results from this being a special case of using a conjugate prior. If you used as a prior random variable with other distribution than Gaussian, there wouldn’t be a simple, closed-form solution for the posterior distribution.

You may want to read more on function convolutions. Random variables are functions.