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

Question

考虑一个非常基本的直线 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 中是否有类似的东西，或者甚至是类似的东西的包？

Answer 1

您可以使用numpy array broadcasting 一步创建您的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 值.