Python实现的威尔逊置信区间算法？

Question

阅读完如何不按平均评分排序后，我想知道是否有Python实现Bernoulli参数的Wilson得分置信区间的下限？

- Jeff Bauer

---

如果 np-cap(1-p-cap) 低于某个阈值，比如30-35，为了更精确，我会使用自由度为(pos+neg)-2的t分布，而不是正态分布。无论如何，这只是我的个人意见。 - luke14free

Answer 1

Reddit使用威尔逊得分区间来对评论进行排名，此处可以找到解释和Python实现链接。

#Rewritten code from /r2/r2/lib/db/_sorts.pyx

from math import sqrt

def confidence(ups, downs):
    n = ups + downs

    if n == 0:
        return 0

    z = 1.0 #1.44 = 85%, 1.96 = 95%
    phat = float(ups) / n
    return ((phat + z*z/(2*n) - z * sqrt((phat*(1-phat)+z*z/(4*n))/n))/(1+z*z/n))

- Steef

---

如果你只是想发布一个链接，请在评论中发布。如果你要将它作为答案发布，提供更多内容或提取代码，这样不需要每个人都跟着链接走，即使链接失效，答案也有价值。 - agf

这个答案应该被修正以包括下面的修改！ - Vladtn

@Vladtn，我刚刚更新了它，使用了Gullevek的答案。如果还有其他问题，请告诉我。 - Amelio Vazquez-Reina

我想补充一下，对于95%的置信区间，z得分应该是1.96而不是1.6。 - Wesley

@Wesley 是的，我认为 1.0 = 85% 也是错误的，已经更新了答案... 这里有一个值表格 http://www.dummies.com/how-to/content/checking-out-statistical-confidence-interval-criti.html - Anentropic

Answer 2

我认为这个答案中的威尔逊调用有误，因为如果你有1个顶和0个踩，你会得到NaN，因为你无法对负值进行sqrt运算。

正确的答案可以在文章如何不按平均分进行排序的Ruby示例中找到：

return ((phat + z*z/(2*n) - z * sqrt((phat*(1-phat)+z*z/(4*n))/n))/(1+z*z/n))

- Gullevek

---

"我认为这个有一个错误的威尔逊调用" - 哪一个？ - onestop

Answer 3

要获取没有连续性修正的Wilson CI，您可以使用statsmodels.stats.proportion中的proportion_confint。 要获取带有连续性修正的Wilson CI，您可以使用以下代码。

# cf. 
# [1] R. G. Newcombe. Two-sided confidence intervals for the single proportion, 1998
# [2] R. G. Newcombe. Interval Estimation for the difference between independent proportions:        comparison of eleven methods, 1998

import numpy as np
from statsmodels.stats.proportion import proportion_confint

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
def propci_wilson_cc(count, nobs, alpha=0.05):
    # get confidence limits for proportion
    # using wilson score method w/ cont correction
    # i.e. Method 4 in Newcombe [1]; 
    # verified via Table 1
    from scipy import stats
    n = nobs
    p = count/n
    q = 1.-p
    z = stats.norm.isf(alpha / 2.)
    z2 = z**2   
    denom = 2*(n+z2)
    num = 2.*n*p+z2-1.-z*np.sqrt(z2-2-1./n+4*p*(n*q+1))    
    ci_l = num/denom
    num = 2.*n*p+z2+1.+z*np.sqrt(z2+2-1./n+4*p*(n*q-1))
    ci_u = num/denom
    if p == 0:
        ci_l = 0.
    elif p == 1:
        ci_u = 1.
    return ci_l, ci_u


def dpropci_wilson_nocc(a,m,b,n,alpha=0.05):
    # get confidence limits for difference in proportions
    #   a/m - b/n
    # using wilson score method WITHOUT cont correction
    # i.e. Method 10 in Newcombe [2]
    # verified via Table II    
    theta = a/m - b/n        
    l1, u1 = proportion_confint(count=a, nobs=m, alpha=0.05, method='wilson')
    l2, u2 = proportion_confint(count=b, nobs=n, alpha=0.05, method='wilson')
    ci_u = theta + np.sqrt((a/m-u1)**2+(b/n-l2)**2)
    ci_l = theta - np.sqrt((a/m-l1)**2+(b/n-u2)**2)     
    return ci_l, ci_u


def dpropci_wilson_cc(a,m,b,n,alpha=0.05):
    # get confidence limits for difference in proportions
    #   a/m - b/n
    # using wilson score method w/ cont correction
    # i.e. Method 11 in Newcombe [2]    
    # verified via Table II  
    theta = a/m - b/n    
    l1, u1 = propci_wilson_cc(count=a, nobs=m, alpha=alpha)
    l2, u2 = propci_wilson_cc(count=b, nobs=n, alpha=alpha)    
    ci_u = theta + np.sqrt((a/m-u1)**2+(b/n-l2)**2)
    ci_l = theta - np.sqrt((a/m-l1)**2+(b/n-u2)**2)     
    return ci_l, ci_u


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
# single proportion testing 
# these come from Newcombe [1] (Table 1)
a_vec = np.array([81, 15, 0, 1])
m_vec = np.array([263, 148, 20, 29])
for (a,m) in zip(a_vec,m_vec):
    l1, u1 = proportion_confint(count=a, nobs=m, alpha=0.05, method='wilson')
    l2, u2 = propci_wilson_cc(count=a, nobs=m, alpha=0.05)
    print(a,m,l1,u1,'   ',l2,u2)

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
# difference in proportions testing 
# these come from Newcombe [2] (Table II)
a_vec = np.array([56,9,6,5,0,0,10,10],dtype=float)
m_vec = np.array([70,10,7,56,10,10,10,10],dtype=float)
b_vec = np.array([48,3,2,0,0,0,0,0],dtype=float)
n_vec = np.array([80,10,7,29,20,10,20,10],dtype=float)

print('\nWilson without CC')
for (a,m,b,n) in zip(a_vec,m_vec,b_vec,n_vec):
    l, u = dpropci_wilson_nocc(a,m,b,n,alpha=0.05)
    print('{:2.0f}/{:2.0f}-{:2.0f}/{:2.0f} ; {:6.4f} ; {:8.4f}, {:8.4f}'.format(a,m,b,n,a/m-b/n,l,u))

print('\nWilson with CC')
for (a,m,b,n) in zip(a_vec,m_vec,b_vec,n_vec):
    l, u = dpropci_wilson_cc(a,m,b,n,alpha=0.05)
    print('{:2.0f}/{:2.0f}-{:2.0f}/{:2.0f} ; {:6.4f} ; {:8.4f}, {:8.4f}'.format(a,m,b,n,a/m-b/n,l,u))

HTH

- batesbatesbates

Answer 4

目前被接受的解决方案似乎使用硬编码的z值（性能最佳）。

如果你想要一个直接的Python等价物来自博客文章中基于置信区间动态计算z值的Ruby公式：

import math

import scipy.stats as st


def ci_lower_bound(pos, n, confidence):
    if n == 0:
        return 0
    z = st.norm.ppf(1 - (1 - confidence) / 2)
    phat = 1.0 * pos / n
    return (phat + z * z / (2 * n) - z * math.sqrt((phat * (1 - phat) + z * z / (4 * n)) / n)) / (1 + z * z / n)

- Loisaida Sam Sandberg

Answer 5

如果您想直接从置信区间计算z并避免安装numpy / scipy，您可以使用以下代码片段：

import math

def binconf(p, n, c=0.95):
  '''
  Calculate binomial confidence interval based on the number of positive and
  negative events observed.  Uses Wilson score and approximations to inverse
  of normal cumulative density function.

  Parameters
  ----------
  p: int
      number of positive events observed
  n: int
      number of negative events observed
  c : optional, [0,1]
      confidence percentage. e.g. 0.95 means 95% confident the probability of
      success lies between the 2 returned values

  Returns
  -------
  theta_low  : float
      lower bound on confidence interval
  theta_high : float
      upper bound on confidence interval
  '''
  p, n = float(p), float(n)
  N    = p + n

  if N == 0.0: return (0.0, 1.0)

  p = p / N
  z = normcdfi(1 - 0.5 * (1-c))

  a1 = 1.0 / (1.0 + z * z / N)
  a2 = p + z * z / (2 * N)
  a3 = z * math.sqrt(p * (1-p) / N + z * z / (4 * N * N))

  return (a1 * (a2 - a3), a1 * (a2 + a3))


def erfi(x):
  """Approximation to inverse error function"""
  a  = 0.147  # MAGIC!!!
  a1 = math.log(1 - x * x)
  a2 = (
    2.0 / (math.pi * a)
    + a1 / 2.0
  )

  return (
    sign(x) *
    math.sqrt( math.sqrt(a2 * a2 - a1 / a) - a2 )
  )


def sign(x):
  if x  < 0: return -1
  if x == 0: return  0
  if x  > 0: return  1


def normcdfi(p, mu=0.0, sigma2=1.0):
  """Inverse CDF of normal distribution"""
  if mu == 0.0 and sigma2 == 1.0:
    return math.sqrt(2) * erfi(2 * p - 1)
  else:
    return mu + math.sqrt(sigma2) * normcdfi(p)

- duckworthd

---

print(binconf(50, 100)) => (0.26291792852889806, 0.41206457669597374) … 50个正事件，100个总事件给出了一个范围，其中上限低于0.5？ - Thomas David Baker

Answer 6

这里是Wilson分数置信区间的简化版本（不需要numpy），并且稍作改进（k值为0和n值时不会导致数学域错误）。该版本带有连续性校正，源代码来自于batesbatesbates在另一个答案中编写的原始代码。还提供了一种纯Python非numpy非连续性校正版本，其中包括两种等效的计算方式（可以通过eqmode参数进行切换，但两种方式都给出了完全相同的非连续性校正结果）。

import math

def propci_wilson_nocc(k, n, z=1.96, eqmode=0):
    # Calculates the Binomial Proportion Confidence Interval using the Wilson Score method without continuation correction
    # Equations eqmode == 1 from: https://en.wikipedia.org/w/index.php?title=Binomial_proportion_confidence_interval&oldid=1101942017#Wilson_score_interval
    # Equations eqmode == 0 from: https://www.evanmiller.org/how-not-to-sort-by-average-rating.html
    # The results should be close to:
    #    from statsmodels.stats.proportion import proportion_confint
    #    proportion_confint(k, n, alpha=0.05, method='wilson')
    #z=1.44 = 85%, 1.96 = 95%
    if n == 0:
        return 0
    p_hat = float(k) / n
    z2 = z**2
    if eqmode == 0:
        ci_l = (p_hat + z2/(2*n) - z*math.sqrt(max(0.0, (p_hat*(1 - p_hat) + z2/(4*n))/n))) / (1 + z2 / n)
    else:
        ci_l = (1.0 / (1.0 + z2/n)) * (p_hat + z2/(2*n)) - (z / (1 + z2/n)) * math.sqrt(max(0.0, (p_hat*(1 - p_hat)/n + z2/(4*(n**2)))))
    if eqmode == 0:
        ci_u = (p_hat + z2/(2*n) + z*math.sqrt(max(0.0, (p_hat*(1 - p_hat) + z2/(4*n))/n))) / (1 + z2 / n)
    else:
        ci_u = (1.0 / (1.0 + z2/n)) * (p_hat + z2/(2*n)) + (z / (1 + z2/n)) * math.sqrt(max(0.0, (p_hat*(1 - p_hat)/n + z2/(4*(n**2)))))
    return [ci_l, ci_u]

def propci_wilson_cc(n, k, z=1.96):
    # Calculates the Binomial Proportion Confidence Interval using the Wilson Score method with continuation correction
    # i.e. Method 4 in Newcombe [1]: R. G. Newcombe. Two-sided confidence intervals for the single proportion, 1998;
    # verified via Table 1
    # originally written by batesbatesbates https://dev59.com/gGkw5IYBdhLWcg3wPYMP#74021634
    p_hat = k/n
    q = 1.0-p
    z2 = z**2
    denom = 2*(n+z2)
    num = 2.0*n*p_hat + z2 - 1.0 - z*math.sqrt(max(0.0, z2 - 2 - 1.0/n + 4*p_hat*(n*q + 1)))
    ci_l = num/denom
    num2 = 2.0*n*p_hat + z2 + 1.0 + z*math.sqrt(max(0.0, z2 + 2 - 1.0/n + 4*p_hat*(n*q - 1)))
    ci_u = num2/denom
    if p_hat == 0:
        ci_l = 0.0
    elif p_hat == 1:
        ci_u = 1.0
    return [ci_l, ci_u]

请注意，返回的值始终会在 [0.0, 1.0] 范围内（因为 p_hat 是 k/n 的比率），这就是为什么它是一个分数而不是真正的置信区间。但是，通过将 ci_l * n 和 ci_u * n 相乘，可以轻松地将其投影回置信区间，这些值将与 k 在相同的域中，并且可以与之并列绘制。

- gaborous

Answer 7

以下是如何计算Wilson Score区间（不需要连续性修正）的更易读版本，由Bartosz Mikulski提供的Python示例：

from math import sqrt
def wilson(p, n, z = 1.96):
    denominator = 1 + z**2/n
    centre_adjusted_probability = p + z*z / (2*n)
    adjusted_standard_deviation = sqrt((p*(1 - p) + z*z / (4*n)) / n)
    
    lower_bound = (centre_adjusted_probability - z*adjusted_standard_deviation) / denominator
    upper_bound = (centre_adjusted_probability + z*adjusted_standard_deviation) / denominator
    return (lower_bound, upper_bound)

- gaborous