2D gaussian distribution does not sum to one? 2D gaussian distribution does not sum to one? numpy numpy

2D gaussian distribution does not sum to one?


NOTE: As stated in the comments bellow, this solution is only valid if you are trying to build a gaussian convolution kernel (or gaussian filter) for image processing purposes. It is not a properly normalized gaussian density function, but it is the form that is used to remove gaussian noise from images.


You are missing the L1 normalization:

ker /= np.abs(ker).sum()

Which will make your kernel behave like an actual density function. Since the grid you have can vary a lot in the magnitude of its values, the above normalization step is needed.

In fact, the gaussian nornalization constant you have could be ommited to just use the L1 norm above. If I'm not worng, you are trying to create a gaussian convolution, and th above is the usual normalization tecnique applied to it.

Your second mistake, as @Praveen has stated, is that you need to sample the grid from [-U//2, U//2]. You can do that as:

i, j = np.mgrid[-U//2:U//2+1, -U//2:U//2+1]

Last, if what you are trying to do is to build a gaussian filter, the size of the kernel is usually estimated from sigma (to avoid zeroes far from the centre) as U//2 <= t * sigma, where t is a truncation parameter usually set t=3 or t=4.


Currently, the (much zoomed-in) contour plot of your ker looks like this:Contour plot of current kernel

As you can see, this looks nothing like a Gaussian kernel. Most of your function dies off from 0 to 1. Looking at the kernel itself reveals that all values do indeed die off really really fast:

>>> ker[0:5, 0:5]array([[  1.592e+001,   3.070e-021,   2.203e-086,   5.879e-195,   0.000e+000],       [  3.070e-021,   5.921e-043,   4.248e-108,   1.134e-216,   0.000e+000],       [  2.203e-086,   4.248e-108,   3.048e-173,   8.136e-282,   0.000e+000],       [  5.879e-195,   1.134e-216,   8.136e-282,   0.000e+000,   0.000e+000],       [  0.000e+000,   0.000e+000,   0.000e+000,   0.000e+000,   0.000e+000]])

The sum value of 15.915 that you get is basically just ker[0, 0]. What all this tells you is that you aren't constructing your grid correctly.

Remember that when creating the kernel on a computer, you'll have to sample it at appropriate points. Sampling it too coarsely will result in your sum not coming out right.

So firstly, if you want a full density centered around mu=0, you'll have to take i and j from -U // 2 to U // 2. But to solve your resolution problem, I recommend taking U number of points between -0.5 and 0.5.

import numpy as npimport matplotlib.pyplot as pltU = 60m = np.linspace(-0.5, 0.5, U)    # 60 points between -1 and 1delta = m[1] - m[0]              # delta^2 is the area of each grid cell(x, y) = np.meshgrid(m, m)       # Create the meshsigma = 0.1norm_constant = 1 / (2 * np.pi * sigma**2)rhs = np.exp(-.5 * (x**2 + y**2) / sigma**2)ker = norm_constant * rhsprint(ker.sum() * delta**2)plt.contour(x, y, ker)plt.axis('equal')plt.show()

This yields a sum close to 1.0, and a kernel centered at mu=0 as expected.Contour plot of corrected kernel

Knowing what range to choose (-0.5 to 0.5) in this case, depends on your function. For instance, if you now take sigma = 2, you'll find that your sum won't work out again, because now you're sampling too finely. Setting your range to be a function of your parameters - something like -5 * sigma to 5 * sigma - might be the best option.