如何在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数组广播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值.