如何在Python中有效地生成具有随机斜率和截距的直线？

Fred S 发表于 Dev

弗雷德·S

考虑一个非常基本的直线蒙特卡罗模拟y = m * x + b，例如为了可视化参数m和不确定度的影响b。m并且b都从正态分布中采样。来自MATLAB背景，我会这样写：

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(start=0, stop=5, step=0.1)

n_data = len(x)
n_rnd = 1000

m = np.random.normal(loc=1, scale=0.3, size=n_rnd) 
b = np.random.normal(loc=5, scale=0.3, size=n_rnd)

y = np.zeros((n_data, n_rnd))  # pre-allocate y

for realization in xrange(n_rnd):
    y[:,realization] = m[realization] * x + b[realization]

plt.plot(x, y, "k", alpha=0.05);

直线的蒙特卡洛模拟

这确实产生了所需的输出，但是我有点觉得必须有一种更“ Pythonic”的方式来做到这一点。我错了吗？如果不是，那么谁能为我提供一些代码示例来说明如何更有效地执行此操作？

举个例子，我正在寻找：在MATLAB中，无需使用循环就可以轻松编写该代码bsxfun()。Python中是否有类似的东西，甚至还有类似的软件包？

物理学家

您可以使用numpy数组广播y一步一步创建数组，如下所示。

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(start=0, stop=5, step=0.1)

n_data = len(x)
n_rnd = 1000

m = np.random.normal(loc=1, scale=0.3, size=n_rnd) 
b = np.random.normal(loc=5, scale=0.3, size=n_rnd)

y = m * x[:, np.newaxis] + b

for val in y.transpose():
    plt.plot(x, val, alpha=0.05)

# Or without the iteration:
# plt.plot(x, y, alpha=0.05)

plt.show()

x[:, np.newaxis]强制x变为形状的列矢量，这与广播(50, 1)相反(50,)。

然后，您可以直接在numpy数组上进行迭代（而不是在其索引上进行迭代），但是必须对数组进行转置（使用y.transpose()），否则对于每次迭代，您将获得每个1000个随机数的x值。

本文收集自互联网，转载请注明来源。

如有侵权，请联系 [email protected] 删除。