Python:计算超几何分布的置信区间

Question

我找到了这个网站:http://www.cluster-text.com/confidence_interval.php计算超几何分布的置信区间，链接到维基百科的 Clopper-Pearson 区间部分:https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval#Clopper-Pearson_interval

我现在想知道如何在 Python 中进行相同的计算。我假设“句子”下方的两个大表达式“Beta 分布又与 F 分布相关，因此可以使用 F 分位数编写 Clopper-Pearson 区间的第三个公式:”是上下边界.出于某种原因，该公式似乎根本不关心人口规模......

维基百科给出的公式调用了一个F函数，我以为是scipy.stats.f中的一些函数。我编写了一个程序来尝试所有这些，我发现使用 stats.f.ppf 作为“神秘的”F 函数可以得到最好的结果。

我设法编写的函数给出的结果与网站给出的结果相当接近，但又不太接近(可能是因为忽略了人口规模)。我需要做什么来解决这个问题并获得正确的置信区间？例如:根据网站，如果人口规模为 80，样本规模为 16，您发现 8 次成功，则 95% 置信区间为 (27.5, 72.5)。但是，我的函数给出 (24.651011149057535, 75.348988850942462)。相当接近，但不完全是。我在代码中包含了另外两个示例。

这是我的代码。我已经包括了(也许是愚蠢的)蛮力尝试来找到要使用的正确 F 函数。 [我不是专业的程序员，也不是专业的统计学家……我尽力了，但很难]

from scipy import stats

def hypergeom_ci(population_size, sample_size, successes, p):
    n = sample_size
    x = successes
    alpha = 1-p
    lower = (1+(n-x+1)/(x*F(alpha/2, 2*x, 2*(n-x+1))))**(-1)
    upper = (1+(n-x)/((x+1)*F(1-alpha/2, 2*(x+1), 2*(n-x))))**(-1)
    return 100 * lower, 100 * upper

population_size, sample_size, successes, p = 80, 16, 8, 0.95

for pk in dir(stats.f): # Try all functions in stats.f and print results
    try:
        F = eval('stats.f.'+pk)
        print(pk, hypergeom_ci(population_size, sample_size, successes, p), sep=' -> ')
    except:
        pass

pk = 'ppf' # I found that stats.f.ppf gives the closest and most reasonable results
F = eval('stats.f.'+pk)

attempts = [(population_size, sample_size, successes, p), (1000000, 70, 8, 0.95), (8883, 70, 8, 0.92)]
h3 = (5.583699200720477, 20.17336485421592)
h1 = (27.5, 72.5)
h2 = (5.0654, 21.2824)

sitegives = [h1, h2, h3] # Answers given by the calculator at http://www.cluster-text.com/confidence_interval.php

for trial, answer in zip(attempts, sitegives):
    print(trial, hypergeom_ci(*trial), answer)

Answer 1

这不是一个直接的答案，但我已经制作了一个计算超几何分布的置信区间的包。还是有点bug，大家使用请慎用。

Github:https://github.com/KeisukeNagakawa/cisim/blob/master/README.rst

PyPI:https://pypi.org/project/cisim/

我也想要您所描述的确切解决方案。