logo标签列表

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

pythonpython-3.xstatistics
3
我发现了这个网站：http://www.cluster-text.com/confidence_interval.php ，它可以计算超几何分布的置信区间，它链接到维基百科的Clopper-Pearson区间部分，链接如下：https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval#Clopper-Pearson_interval 我现在想知道如何在Python中进行同样的计算。我假设在“The beta distribution is, in turn, related to the F-distribution so a third formulation of the Clopper-Pearson interval can be written using F quantiles:” 下面的两个大表达式是下限和上限。由于某些原因，该公式似乎完全不关心人口数量...
维基百科中介绍的公式调用了一个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)
1
这个差异不就是因为维基百科页面上的Clopper-Pearson区间是针对二项式分布(有放回抽样)而不是超几何分布(无放回抽样)吗？这也解释了为什么人口数量没有进入计算。此外，查看该网站的Javascript源代码，它似乎与维基百科部分完全无关。 - Stefan Zobel
1个回答
0
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接