Numpy.random.normal会给出不好的结果

Question

我试图用numpy.random.normal来模拟一个随机数。根据该随机数（平均值=0，标准值=1）

1. 我绘制了多个大小相似的样本（例如，m=100）
2. 我计算每个样本的标准偏差
3. 我取所有标准偏差的平均值

理论统计，还有R告诉我，这必须收敛到所选的std（即1）。但不知何故，使用numpy（和scipy.stats）时，它并没有这样做

此代码生成一个显示此奇怪行为的图形：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, tstd

# system setup
m = 100         # number of measurments
sigma = 1       # sensor std

ez = np.arange(1,6,.05)
sample_sizes = [int(10**e) for e in ez]

# testing normal and std - they seem to work fine
sig_est = []
for n in sample_sizes:
    sample = np.random.normal(0, sigma, (n*m))
    sig_est += [np.std(sample)]
plt.plot(ez, sig_est, marker='.', color='b', ls='', label='numpy - no means')

# numpy implementation of problem
sig_est = []
for n in sample_sizes:
    sample = np.random.normal(0, sigma, (n,m))
    sigma_est = np.std(sample, axis=1)
    sig_est += [np.mean(sigma_est)]
plt.plot(ez, sig_est, marker='.', color='k', ls='', label='numpy')

# scipy.stats implementation
sig_est = []
for n in sample_sizes:
    sample = norm.rvs(loc=0, scale=sigma, size=(n,m))
    sigma_est = tstd(sample, axis=1)
    sig_est += [np.mean(sigma_est)]
plt.plot(ez, sig_est, marker='.', color='r', ls='', label='scipy.stats')

plt.gca().set(xlabel = 'Number of samples [log10]')
plt.gca().legend()
plt.gca().grid(color='.9')
plt.show()

有什么想法吗

Answer 1

这是一个有趣的问题，因为这不是一个随机数生成器问题，而是一个数学问题：-）简单的答案是，一切都按预期进行

重点是，在第一个示例中，您正在获取越来越大的i.i.d.高斯样本，并使用np.std计算其标准偏差。这会收敛到1，如图所示

在第二个图中，您计算的标准偏差始终超过100个元素，然后对这些元素求平均值。通过这种方式，您不是计算许多元素的极限std，而是计算标准偏差估计值的偏差。正如你发现的，这不是零！这有两个原因：

• 标准偏差的默认numpy实现是方差估计量的平方根，该方差最小化二次风险（即二次误差的1/n和）。这是而不是方差的无偏估计量，它将从1/（n-1）开始。您可以通过将参数ddof=1传递给np.std来获得后者，请参见此处的文档：https://numpy.org/doc/stable/reference/generated/numpy.std.html
• 。。。但即使你这样做了，你也不会得到零偏差。这是因为你绘制的是std，而不是方差；i、 为了得到精确的1，你应该在计算np.std和取平均值之前将结果平方。你可以看到，如果你更换你的线路

sig_est += [np.mean(sigma_est)]  # equivalent to sig_est.append(np.mean(sigma_est))

借

sig_est.append(np.mean(np.std(sample, axis=1, ddof=1)**2))

在代码的第二个块中，您确实会收敛到1

至于最后一个使用scipy的实现，它似乎使用了另一个规范化：https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.tstd.html

他们称之为“无偏”，但它显然不是，一方面是因为你的图清楚地显示了它，另一方面是因为获得无偏估计量（高斯）的确切因子比n/（n-1）复杂得多，请参见这里：https://en.wikipedia.org/wiki/Unbiased_estimation_of_standard_deviation