Python 计算置信区间
用Python实现http://vassarstats.net/prop1.html 所计算的置信区间
import math
#计算置信区间的函数
def calc(r,n):
if n < r :
print ('r cannot be greater than n.')
return
if math.floor(r) < r :
print ('r must be an integer value.')
return
if math.floor(n) < n :
print ('n must be an integer value.')
return
p = round((r/n)*10000)/10000
print ('p',p)
q = 1-p
p = round(p*10000)/10000
print ('p',p)
z=1.95996
zsq = z*z
# l95a
num = (2*n*p)+zsq-(z*math.sqrt(zsq+(4*n*p*q)))
denom = 2*(n+zsq)
l95a = num/denom
if p == 0 :
l95a = 0
l95a = round(l95a*10000)/10000
# u95a
num = (2*n*p)+zsq+(z*math.sqrt(zsq+(4*n*p*q)))
denom = 2*(n+zsq)
u95a = num/denom
print ('u95a',u95a)
if p == 1 :
u95a = 1
u95a = round(u95a*10000)/10000
print('no continuity correction',l95a,'-',u95a)
# l95b
num = (2*n*p)+zsq-1-(z*math.sqrt(zsq-2-(1/n)+4*p*((n*q)+1)))
denom = 2*(n+zsq)
l95b = num/denom
if p==0 :
l95b = 0
l95b = round(l95b*100000000)/100000000
# u95b
num = (2*n*p)+zsq+1+(z*math.sqrt(zsq+2-(1/n)+4*p*((n*q)-1)))
denom = 2*(n+zsq)
u95b = num/denom
if p==1 :
u95b = 1
u95b = round(u95b*100000000)/100000000
print('including continuity correction',l95b,'-',u95b)
return
#测试k = 12 n = 40
calc(12,40)
print('THE END')