binomial_confidence_plot.py

# http://slobin.livejournal.com/515819.html 

import math

def interval(k, n):
    # http://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval
    # Wilson score interval

    z = 1.96 #95% confidence
    p = 1.0*k / n
    a = p + (z*z)/(2.0*n)
    b = z*(math.sqrt(p*(1-p)/n  + z*z/(4*n*n)))
    c = 1/(1.0 + z*z/n)
    return (a - b)*c, (a + b)*c



def main():
    #print interval(10, 20)
    data = [
        (2003,   1),
        (2004,   6),
        (2005,  18),
        (2006,  92),
        (2007, 269),
        (2008, 257),
        (2009, 285),
        (2010, 341),
        (2011, 296),
        (2012, 308),
        (2013, 318),
        (2014, 499),
        (2015, 286),
    ]

    n = 2976
    for year, k in data:
        a, b = interval(k, n)
        p = 1.0*k/n
        assert a < p < b
        ia = int(a * 300.0)
        ib = int(b * 300.0)
        ip = int(p * 300.0)
        ib = max(ib, ia+1) # do not put [] to one point
        assert ia <= ip <= ib, (ia, ip, ib)
        line = [' '] * (ib + 1)
        line[0:ip] = '*'*ip
        line[ib] = ']'
        line[ia] = '['
        print '%s % 4d %s' % (year, k, ''.join(line))


main()