2012-04-05 2 views
24

How Not to Sort by Average Rating을 읽은 후 누구나 Bernoulli 매개 변수에 대한 Wilson 점수 신뢰 구간의 하위 경계를 Python으로 구현했는지 궁금합니다.파이썬에서 Wilson 점수 간격을 구현합니까?

+0

보다 정밀도 경우 N 개의 * P는 캡 * (1-p-cap)이 일정한 임계 값, 예를 들어 30-35보다 작 으면 정상 분배기 대신 df : (pos + neg) -2를 사용하는 t- 분배를 사용합니다. 아무리 해도. 그래도 내 두 센트 – luke14free

답변

26

레딧 코멘트는 설명과 파이썬 구현은 당신이 아래로 당신이 얻을 1 0이있는 경우 때문에 나는이 하나가 잘못 윌슨 전화를 가지고 있다고 생각

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

당신은 단지 링크를 게시하려고하는 경우, 코멘트에서 해. 답으로 게시하는 경우 콘텐츠에서 더 많은 정보를 제공하거나 코드를 추출하여 모든 사람이 링크를 따라야하는 것은 아니므로 링크가 작동하지 않아도 대답은 가치가 있습니다. – agf

+0

그의 답변을 받아 들일 수 있도록 Steef의 링크에서 코드를 추가했습니다. –

+2

이 답변은 아래 수정을 포함하도록 수정해야합니다! – Vladtn

18

here을 찾을 수 있습니다 순위에 대한 윌슨 점수 간격을 사용 NaN 음수 값으로 sqrt을 수행 할 수 없기 때문입니다.

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

실제로 바인딩 신뢰에서 직접 z를 계산하기 좋아하고 설치하지 않도록 할 것입니다 경우 NumPy와/scipy : 기사 How not to sort by average page에서 루비 예를 볼 때

올바른

찾을 수 있습니다 , 당신은 당신이 statsmodels.stats.proportionproportion_confint을 사용할 수 있습니다, 연속성 보정없이 윌슨 CI를 얻으려면,

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

을 코드의 다음 코드를 사용할 수 있습니다. 연속성 보정 기능이있는 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

2

허용되는 용액 (최고 성능) 하드 코딩 된 Z 값을 사용하는 것으로 보인다. 동적 Z 값과 the blogpost에서 루비 화학식 직접 파이썬 상당 원하는 경우에

는 (신뢰 구간 기준) :

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)