从Matlab到Python的算法等价性

Question

我用下面这个小的m-file在Matlab中绘制了一个3D网格：

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

我将通过使用Python及其相应的模块，使用以下代码段来获得相同的结果：

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()

但是Python的输出与Matlab的输出相比有很大的不同，事实上是无法接受的。我担心会出现函数误用，就像linespace、enumerate或mgrid一样...

有没有人有什么想法？...

PS：不幸的是，我不能在这个帖子中插入输出图形...！

最好

- user1393214

---

tt = np.linspace(0, 19, 20) 看起来应该是201而不是19吗？ - Brian L

我根据您的建议更新了问题... - user1393214

Answer 1

根据我的理解，等效的解决方案应该是：

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()

不幸的是，使用matplotlib进行3D绘图并不像在matlab中那样直截了当。

以下是来自此代码的图表：

- swenzel

---

感谢您的贡献...但是有两个问题： 1- 结果与我的解决方案和Matlab驱动代码不同。您如何评估？ 2- 我不明白您如何合法地将“sigma”替换为“np.sqrt(2)*n”，因为它应该是1...您能否请给予一些评论？ - user1393214

@Ordenador 如果你忽略我发布的图表并使用我提供的代码，它将产生与你的Matlab脚本完全相同的结果，即 sigma = sqrt(2)./n ;) 然而，你的Python脚本对 norm(t*n, 0, 1) 在t = 0到x上进行积分，这与 normcdf(x, 0, n) 相同。因此，你的Python脚本中的sigma实际上是 n 而不是 1，因此你的脚本在sigma的值上并不一致，我不知道你实际需要什么。我选择将 sigma = np.sqrt(2)*n 作为我的图表，因为出于某种原因我认为这是一个很好的折衷方案... - swenzel

抱歉，我刚意识到我搞混了一些东西... 对于t从0到x的norm(t*n,0,1)积分等同于normcdf(t,0,1/n)，它看起来与normcdf(t,0,sqrt(2)/n)相似，这就是为什么你的脚本看起来产生了相同的图形，但从技术上讲它并不是。我更新了图形以显示我的代码生成的内容。 - swenzel

Answer 2

您的Matlab代码通过x生成了201个点：

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

当你的Python代码只生成20个点时：

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

也许它会导致准确性问题？ 尝试使用以下代码：

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

- akarilimano

Answer 3

解决问题的关键点包括：

1- 对于提供的linespace和mgrid函数相关维度的等效性的必要性...

2- 利用更高密度的网格使路径尽可能平滑...

3- 应用三维绘图函数，例如plot_surf...

当前代码完全有效...

- user1393214