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

Question

我需要从方差为 1 和我指定的真实 mu(平均值)的正态分布中生成大小为 200 (n=200) 的随机样本；然后，我根据假设测试抽签结果:mu <= 1。我需要对 400 个潜在的真实 theta 中的每一个执行此操作，并且对于每个真实的 theta，我需要将其复制 200 次。

我已经对 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 = 1600 万 个随机值(消耗约 122 MB 内存；使用 draws.nbytes/1024/1024)，并使用 SciPy 对每个 θ 值的 200 个观测值的 200 个样本中的每一个样本运行单侧单样本 t 检验:

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$')