bivariate normal distribution function buggy?

Community matplotlib-users
1

Hi,

I think the bivariate normal probability distribution function implementation in Matplotlib is buggy. See the example code below, and note that the change in the x directions and the y directions is not uniform as it should be.

The implementation (in mlab.py) is appended.

                                                                Faheem.

···

***************************************************************************
minx, maxx = 0, 13
miny, maxy = 0, 13
X, Y = pylab.meshgrid(range(minx, maxx), range(miny, maxy))
conc = 1000*pylab.bivariate_normal(X, Y, sigmax=5.0, sigmay=5.0, sigmaxy=0, mux=6, muy=6)

[[ 0.08 0.11 0.13 0.15 0.16 0.17 0.17 0.17 0.16 0.15 0.13 0.11 0.08]
[ 0.25 0.32 0.38 0.44 0.48 0.51 0.52 0.51 0.48 0.44 0.38 0.32 0.25]
[ 0.63 0.78 0.93 1.07 1.19 1.26 1.29 1.26 1.19 1.07 0.93 0.78 0.63]
[ 1.26 1.57 1.88 2.16 2.39 2.54 2.59 2.54 2.39 2.16 1.88 1.57 1.26]
[ 2.08 2.59 3.1 3.56 3.94 4.18 4.27 4.18 3.94 3.56 3.1 2.59 2.08]
[ 2.8 3.49 4.18 4.81 5.32 5.65 5.76 5.65 5.32 4.81 4.18 3.49 2.8 ]
[ 3.1 3.86 4.62 5.32 5.88 6.24 6.37 6.24 5.88 5.32 4.62 3.86 3.1 ]
[ 2.8 3.49 4.18 4.81 5.32 5.65 5.76 5.65 5.32 4.81 4.18 3.49 2.8 ]
[ 2.08 2.59 3.1 3.56 3.94 4.18 4.27 4.18 3.94 3.56 3.1 2.59 2.08]
[ 1.26 1.57 1.88 2.16 2.39 2.54 2.59 2.54 2.39 2.16 1.88 1.57 1.26]
[ 0.63 0.78 0.93 1.07 1.19 1.26 1.29 1.26 1.19 1.07 0.93 0.78 0.63]
[ 0.25 0.32 0.38 0.44 0.48 0.51 0.52 0.51 0.48 0.44 0.38 0.32 0.25]
[ 0.08 0.11 0.13 0.15 0.16 0.17 0.17 0.17 0.16 0.15 0.13 0.11 0.08]]

def bivariate_normal(X, Y, sigmax=1.0, sigmay=1.0,
                      mux=0.0, muy=0.0, sigmaxy=0.0):
     """
     Bivariate gaussan distribution for equal shape X, Y

     http://mathworld.wolfram.com/BivariateNormalDistribution.html
     """
     Xmu = X-mux
     Ymu = Y-muy

     rho = sigmaxy/(sigmax*sigmay)
     z = Xmu**2/sigmax**2 + Ymu**2/sigmay - 2*rho*Xmu*Ymu/(sigmax*sigmay)
     return 1.0/(2*pi*sigmax*sigmay*(1-rho**2)) * exp( -z/(2*(1-rho**2)))

2

Hello,

It looks like there are a couple of typos in the equations used in the
function.

I think the bivariate normal probability distribution function
implementation in Matplotlib is buggy. See the example code below, and
note that the change in the x directions and the y directions is not
uniform as it should be.
The implementation (in mlab.py) is appended.

[snip]

 z = Xmu\*\*2/sigmax\*\*2 \+ Ymu\*\*2/sigmay \- 2\*rho\*Xmu\*Ymu/\(sigmax\*sigmay\)

                                           ^ **2

 return 1\.0/\(2\*pi\*sigmax\*sigmay\*\(1\-rho\*\*2\)\) \* exp\( \-z/\(2\*\(1\-rho\*\*2\)\)\)

                                      ^ sqrt()

Compare to equations (1) and (2) from the link in the docstring:

I've attached a patch.

Mike

mlab.py.diff (579 Bytes)

···

On Friday 07 July 2006 12:17, Faheem Mitha wrote: