Algorithm equalivence from Matlab to Python

Question

I've plotted a 3-d mesh in Matlab by below little m-file:

[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);

I am going to acquire the same result by utilization of Python and its corresponding modules, by below code snippet:

import numpy as np
from scipy.integrate import quad
import matplotlib.pyplot as plt

sigma = 1

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

tt = np.linspace(0, 20, 2000)
nn = np.linspace(1, 100, 100)  

T = np.zeros([len(tt), len(nn)])

for i,t in enumerate(tt):
    for j,n in enumerate(nn):
        T[i, j], _ = quad(integrand, -np.inf, t, args=(n,))

x, y = np.mgrid[0:20:0.01, 1:101:1]

plt.pcolormesh(x, y, T)

plt.show()

But the output of the Python is is considerably different with the Matlab one, and as a matter of fact is unacceptable. I am afraid of wrong utilization of the functions just like linespace, enumerate or mgrid...

Does somebody have any idea about?!...

PS. Unfortunately, I couldn't insert the output plots within this thread...!

Best

..............................

Edit: I changed the linespace and mgrid intervals and replaced plot_surface method... The output is 3d now with the suitable accuracy and smoothness...

Answer 1

From what I see the equivalent solution would be:

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

x, n = np.mgrid[0:20:0.01, 1:100:1]

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

f = norm.cdf(x, mu, sigma)

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(x, n, f, rstride=x.shape[0]//20, cstride=x.shape[1]//20, alpha=0.3)

plt.show()

Unfortunately 3D plotting with matplotlib is not as straight forward as with matlab.

Here is the plot from this code:

Answer 2

Your Matlab code generate 201 points through x:

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

While your Python code generate only 20 points:

tt = np.linspace(0, 19, 20)

Maybe it's causing accuracy problems? Try this code:

tt = np.linspace(0, 20, 201)