如何在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值。