# Pointwise multiplication of two Gaussians and normalization in Python

The Bayesian update contains a multiplication of the Prior and the Likelihood. The area under the bell curve of a gaussian sums up to 1. We know that this only holds if we integrate from -Inf to +Inf.

This is not feasible on a computer, so I tried to somehow approximate the normalization step after a multiplication of two gaussians.

First, I create one gaussian using the quasi-continuous space between -50.0 and +50.0:

from scipy.stats import norm
import numpy as np

x = np.linspace(-50.0, 50.0, 10000)
y1 = norm.pdf(x)


When I now sum over all values of y1, I get 99.99 - as the sum does not take into account that I have a point each 1/100 of a full step.

> print y1.sum()
> 99.99


When I'd do a Riemann-like integration, like multiplying each x distance with the y value, then it will some up to appx. one:

> print np.multiply(y1, 0.01).sum()
> 0.9999


Now I create another gaussian, slightly shifted:

y2 = np.norm(x+2)


When I plot these two curves, I see the desired and well-known picture:

import matplotlib.pyplot as plt
%matplotlib inline

plt.ylim(0.0,0.5)
plt.xlim(-6.0,6.0)
plt.plot(x, y1)
plt.plot(x, y2)


So, now let's assume that I multiply these two gaussians:

y3 = (y1 * y2)


When I would add this y3 to the plot, of course the values would be too small as the y3 are not normalized yet - but, I also cannot normalize the values to sum up to 1 - as seen above, my original PDF-values of y1 and y2 have either summed up to 99.99 due to the sampling. So, I do an approximation:

y3 = ((y1.sum() + y2.sum()) / 2.0) * (y3.astype(float) / y3.astype(float).sum())


which is kind of a scaled normalization. When I now plot the y3 along with the other two gaussians, I get an image that does not look to be totally wrong:

plt.ylim(0.0,.7)
plt.xlim(-8.0,6.0)
plt.plot(x, y1)
plt.plot(x, y2)
plt.plot(x, y3, '--')


Now, my questions are:

• Is there any other (more pythonic way) to multiply and normalize two gaussians? Or more generally, two unknown PDFs?
• Is the approach that I have implemented methodically o.k.? If not, what are the weaknesses?
• Your question indeed seems to be closer to the content at Cross Validated (CV) than to one at this site. Do you want me to migrate the question to CV? – Aleksandr Blekh Jun 29 '16 at 2:54
• Let me try... You are welcome. – Aleksandr Blekh Jun 29 '16 at 9:50
• If you would just do the arithmetic, you will discover that the product of two Gaussians is another Gaussian. That will save you a great deal of coding and result in a program that is far faster and more accurate. Your other questions are quite broad because they beg for a treatise on numerical integration. – whuber Jun 29 '16 at 12:44
• This is true for gaussians, but I am aiming towards a solution that can handle probability distributions (approximated by the probability mass functions) which are unknown in the first place. You are right, in the case of two gaussians I could easily calculate the new mean and variance for the product of two gaussians using a formulae as stated in johndcook.com/blog/2012/10/29/product-of-normal-pdfs or more precisely in blog.jafma.net/2010/11/09/… – Regenschein Jun 29 '16 at 13:15
• math.stackexchange.com/questions/101062/… – Kevin May 20 '18 at 3:50