如何使用numpy正态分布采样器无循环高效生成M×N正态样本数组

需求为从M个各自带有对应标准差的源值生成服从正态分布的样本，每个源值对应生成N个样本，期望输出形状为M×N的数组。
已知Python层循环可实现该需求，但需要依托numpy原生向量化能力的高效无循环实现方案。

原错误实现 #

import numpy as np

# 初始数据
M = 100
x = np.arange(M)
y = x**2
y_err = y * 0.1

# 每个数据点采样N次
N = 1000
N_samples = np.random.normal(loc=y, scale=y_err, size=N)

运行后触发广播错误：

ValueError: shape mismatch: objects cannot be broadcast to a single shape

报错原因：loc、scale参数均为形状(M,)的一维数组，指定输出size为N时，参数形状无法和目标输出形状对齐，无法完成广播。

numpy的随机数生成函数原生支持广播规则，只需要给均值、标准差参数新增一个长度为1的维度对应采样次数轴，numpy会自动在底层完成广播计算，无需编写Python层循环，运行效率和原生向量化运算一致。

import numpy as np

M = 100
x = np.arange(M)
y = x**2
y_err = y * 0.1
N = 1000

# 为参数新增长度为1的轴，自动广播生成M×N的样本矩阵
samples = np.random.normal(
    loc=y[:, None],  # 形状从(M,)转换为(M, 1)，np.newaxis和None写法等价
    scale=y_err[:, None],
    size=(M, N)
)

结果验证 #

• 执行samples.shape可得到输出(100, 1000)，完全匹配M×N的形状要求
• 逐行计算样本均值、标准差，会分别趋近于对应位置的源值y、标准差y_err，符合正态分布参数预期

内容的提问来源于stack exchange，提问作者ascendants