针对零假设测试 80,000 多个模拟正态分布观察集

Question

我需要从方差为 1 和我指定的真实 mu（平均值）的正态分布中生成大小为 200（n=200）的随机样本；然后，我根据一个假设检验平局：mu

我已经为 n=1 做了这个，但我意识到我的方法是不可复制的。对于每 400 个 theta，我运行以下命令：

sample_r200n1_t2=normal(loc=-0.99, scale=1, size=200)
sample_r200n1_t3=normal(loc=-0.98, scale=1, size=200)
sample_r200n1_t4=normal(loc=-0.97, scale=1, size=200)
sample_r200n1_t5=normal(loc=-0.96, scale=1, size=200)
... on and on to loc = 3 

然后，我分别测试了生成数组中的每个元素。然而，这种方法需要我生成数以万计的样本，我生成与每个样本相关的平均值，然后根据我的标准测试该平均值。这必须完成 80,000 次（除此之外，我还需要针对多个不同大小的 n 执行此操作）。显然 - 这不是要采取的方法。

我怎样才能达到我想要的结果？例如，有没有办法生成一组样本均值并将这些均值放入一个数组中，每个 theta 一个？然后我可以像以前一样测试。或者，还有其他方法吗？

Answer 1

您可以在一个 numpy 数组中生成所有 200*200*400 = 16 million 随机值（这会消耗大约 122 兆字节的内存；请使用 draws.nbytes/1024/1024 进行检查），并使用 SciPy 对每个随机值运行单边单样本 t 检验对于每个 theta 值的 200 个观测值的 200 个样本：

from numpy.random import normal
from scipy.stats import ttest_1samp
import matplotlib.pyplot as plt

# Array of loc values; for each loc, we draw 200 
# samples of 200 normally distributed observations
locs = np.linspace(-1, 3, 401)

# Array of shape (401, 200, 200) = (locs, samples, observations)
# Note that 200 draws of 200 i.i.d. observations is the same as
# 1 draw of 200*200 i.i.d. observations, reshaped to (200, 200)
draws = np.array([normal(loc=x, scale=1, size=200*200)
                  for x in locs]).reshape(401, 200, 200)

# axis=1 computes t-test across columns.
# Alternative hypothesis that sample mean
# is less than the population mean of 1 implies a null
# hypothesis that sample mean is greater than or equal to
# the population mean
tstats, pvals = ttest_1samp(draws, 1, alternative='less', axis=1)

# Count how many out of 200 t-tests reject the null hypothesis
# at the alpha=0.05 level
rejects = (pvals < 0.05).sum(axis=1)

# Visual check: p-values should be low for sample means
# far below 1, as these tests should reject the null 
# hypothesis that sample mean >= 1
plt.plot(locs, rejects)
plt.axvline(1, c='r')
plt.title('Number of t-tests rejecting $H_0 : \mu \geq 1$ with $p < 0.05$')
plt.xlabel('Known sample mean $\\theta$')