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.plot(x, y1)
plt.plot(x, y2)

Two bell curves

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.plot(x, y1)
plt.plot(x, y2)
plt.plot(x, y3, '--')

multiplied gaussians

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?

migrated from datascience.stackexchange.com Jun 29 '16 at 9:50

This question came from our site for Data science professionals, Machine Learning specialists, and those interested in learning more about the field.

  • $\begingroup$ 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? $\endgroup$ – Aleksandr Blekh Jun 29 '16 at 2:54
  • $\begingroup$ Let me try... You are welcome. $\endgroup$ – Aleksandr Blekh Jun 29 '16 at 9:50
  • 1
    $\begingroup$ 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. $\endgroup$ – whuber Jun 29 '16 at 12:44
  • $\begingroup$ 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/… $\endgroup$ – Regenschein Jun 29 '16 at 13:15
  • 1
    $\begingroup$ math.stackexchange.com/questions/101062/… $\endgroup$ – Kevin May 20 '18 at 3:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.