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

Test 80,000+ simulated normal distribution observation sets against a null hypothesis

我需要从方差为 1 和我指定的真实 mu(平均值)的正态分布中生成一个大小为 200 (n=200) 的 运行dom 样本;然后,我根据假设测试平局: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 一个?然后我可以像以前一样测试。或者,还有其他方法吗?

您可以在一个 numpy 数组中生成所有 200*200*400 = 16 million 个随机值(这会消耗大约 122 兆字节的内存;请查看 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$')