0

I have a piece of code which is supposed to plot a Gaussian distribution over bivariate data, using arrays and the imshow() function from pylab. My code looks like:

from pylab import *
x = zeros((101,101))
mean_i = 50
mean_j = 50
sigma = 10
max_ = 100
for i in range(101):
    for j in range(101):
        x[i,j] = max_*exp(((-(i-mean_i)**2)-(j-mean_j)**2)/(2*(sigma**2)))
for i in range(101):
    plot(i,x[i,50], 'rx')
show()
imshow(x)
show()

where my equation for x[i,j] is the product of the Gaussian for i and j.

If you run the code, you will see that two plots are produced; one is a colourmap of the i/j array, which should have a bright spot in the centre and decay away, as a Gaussian would, and the other is a plot of the i values for j=50, so should be a regular Gaussian curve.

This is not what happens. When I run the code, there are steps, like the data is not continuous when I feel like it should be. Why is this? Is it a problem with my equation, or my plotting technique?

Dan Pollard
  • 51
  • 1
  • 7

1 Answers1

1

You are obviously using python 2 and run into the usual problem of integer division

The easiest way out is to make sigma a float, 

sigma = 10.0

You can also cast the denominator to float

x[i,j] = max_*exp(((-(i-mean_i)**2)-(j-mean_j)**2)/float(2*(sigma**2)))

Or you can, at the beginning of your script, import the python 3 division behaviour:

from __future__ import division

and keep the rest the same.

I would also suggest to adapt your code to get rid of the for-loops. Also, importing everything from pylab is not recommended. 

import numpy as np
import matplotlib.pyplot as plt

mean_i = 50
mean_j = 50
sigma = 10.
max_ = 100
X,Y = np.meshgrid(np.arange(101), np.arange(101))
x = max_*np.exp(((-(X-mean_i)**2)-(Y-mean_j)**2)/(2*(sigma**2)))

plt.plot(np.arange(101),x[:,50], 'rx')
plt.show()
plt.imshow(x)
plt.show()
ImportanceOfBeingErnest
  • 321,279
  • 53
  • 665
  • 712