Implementación de Python del intervalo de puntuación de Wilson?

Question

Después de leer How Not to Sort by Average Rating, tenía curiosidad si alguien tiene una implementación en Python de un límite inferior del intervalo de confianza de puntuación de Wilson para un parámetro de Bernoulli?

Answer 1

Reddit utiliza el intervalo de puntuación de Wilson para la clasificación comentario, una implementación de la explicación y la pitón se puede encontrar here

#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)) 

Answer 2

Creo que éste tiene una llamada Wilson mal, porque si usted tiene 1 a 0 abajo se obtiene NaN porque no puede hacer un sqrt en el valor negativo.

el correcto se puede encontrar cuando se mira en el ejemplo del artículo de rubí How not to sort by average page:

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

Answer 3

Si desea calcular z realidad directamente de un límite de confianza y desea evitar la instalación de numpy/scipy , puede utilizar el siguiente fragmento de código,

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) 

Answer 4

Para obtener el CI Wilson y sin corrección de continuidad, se puede utilizar en proportion_confintstatsmodels.stats.proportion. Para obtener el Wilson CI con corrección de continuidad, puede usar el siguiente código.

# 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

Answer 5

La solución aceptada parece usar un valor z-codificado (mejor para el rendimiento).

En el caso de que quería un equivalente pitón directa de la fórmula rubí de the blogpost con un valor z dinámico (basado en el intervalo de confianza):

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)