2

Here is my first steps within the NumPy world. As a matter of fact the target is plotting below 2-D function as a 3-D mesh:

N = \frac{n}{2\sigma\sqrt{\pi}}\exp^{-\frac{n^{2}x^{2}}{4\sigma^{2}}}

That could been done as a piece a cake in Matlab with below snippet:

[x,n] = meshgrid(0:0.1:20, 1:1:100);

mu = 0;
sigma = sqrt(2)./n;

f = normcdf(x,mu,sigma);
mesh(x,n,f);

But the bloody result is ugly enough to drive me trying Python capabilities to generate scientific plots.

I searched something and found that the primary steps to hit above mark in Pyhton might be acquired by below snippet:

from matplotlib.patches import Polygon
import numpy as np
from scipy.integrate import quad
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

sigma = 1

def integrand(x,n):
    return (n/(2*sigma*np.sqrt(np.pi)))*np.exp(-(n**2*x**2)/(4*sigma**2))

t = np.linespace(0, 20, 0.01)
n = np.linespace(1, 100, 1)

lower_bound = -100000000000000000000 #-inf
upper_bound = t

tt, nn = np.meshgrid(t,n)

real_integral = quad(integrand(tt,nn), lower_bound, upper_bound)

Axes3D.plot_trisurf(real_integral, tt,nn)

Edit: With due attention to more investigations on Greg's advices, above code is the most updated snippet.

Here is the generated exception:

RuntimeError: infinity comparisons don't work for you

It is seemingly referring to the quad call...

Would you please helping me to handle this integrating-plotting problem?!...

Best

User
  • 952
  • 2
  • 21
  • 43
  • There is a matplotlib plotting package which is somewhat equal to matplot. It has simple examples on how to plot data. Google for it and plot yourself. Also post if you have tried plotting using any of the other packages. – Versatile May 20 '15 at 19:44
  • @Versatile: Thank you to mention that, but as I said I've checked the samples but most of them are around the simple cases, not 2-D integral functions... – User May 20 '15 at 19:49

1 Answers1

0

Just a few hints to get you in the right direction. numpy.meshgrid can do the same as MatLABs function: http://docs.scipy.org/doc/numpy/reference/generated/numpy.meshgrid.html

When you have x and n you can do math just like in matlab:

sigma = numpy.sqrt(2)/n

(in python multiplication/division is default index by index - no dot needed)

scipy has a lot more advanced functions, see for example How to calculate cumulative normal distribution in Python for a 1D case.

For plotting you can use matplotlibs pcolormesh:

import matplotlib.pyplot as plt
plt.pcolormesh(x,n,real_integral)

Hope this helps until someone can give you a more detailed answer.

Community
  • 1
  • 1
Greg
  • 151
  • 7
  • You are welcome - I just added a bit more information to help you with plotting. – Greg May 20 '15 at 20:57
  • Nice suggestion to use `pcolormesh`... But actually the current problem is just before plotting, i.e. within the integral computation, as I've illustrated within the question... – User May 20 '15 at 20:59
  • A totally different direction, but if you are trying to plot the cumulative normal distribution wouldn't it be easier to use the error function, scipy.special.erf? Alternatively if you want to use numerical integration, numpy.cumsum might be relevant for you. – Greg May 20 '15 at 21:07