Binomial distribution

From Chapter 1 of Information Theory, Inference and Learning Algorithms

Question

A coin has probability f of coming up heads, tossed N times. Probability of getting r heads? Mean and variance of r?

$$P(r|f,N) = {N \choose r} f^r(1-f)^{N-r}$$

Factorial

$${N \choose r} = \frac{N!}{(N-r)!r!}$$

Stirling's approximation for the factorial function $$ x! \simeq x^x e^{-x} $$


In [3]:
def fact(n):
    "This function computes a factorial" 
    if (n <= 1): return 1
    else: return n * fact(n - 1)

def stirling(n):
    "Approximates factorial with Stirling's approximation"
    return( (n**n) * (math.e ** (-n)) )
    
def bin(N=10, r=3, f=0.5):
    Nr = (fact(N)/(fact(N-r)*fact(r)))
    prob = Nr * (f**r) * ((1-f)**(N-r))
    return(prob)

def range2(start, stop, step):
    c = start
    while c < stop:
        c += step
        yield c

def plotfn(fn, start=0, end=100, step=1, ylab="", xlab="Range", titletxt=""):
    rng = lambda: range2(start, end, step)
    plot(list(rng()), [fn(i) for i in rng()])
    ylabel(ylab)
    xlabel(xlab)
    title(titletxt)

Varying across number of coins thrown


In [116]:
plotfn((lambda x: bin(r=x)), start=0, end=10, 
       xlab="Number of heads desired", ylab="Probability", titletxt="Probability of x heads given 10 coins thrown")


Varying probability of one coin being head


In [115]:
plotfn((lambda x: bin(f=x)), start=0, end=1, step=0.05, 
       xlab="Probability of one coin being head", ylab="Probability", titletxt="Probability of 5 heads given 10 throws")


What if I want to calculate the probability of getting five heads or more?


In [121]:
def cumbin(N, r, f):
    return(sum([bin(N, x, f) for x in range(r, N)]))

In [123]:
plotfn((lambda x: cumbin(N=10, r=x, f=0.5)), start=0, end=10, 
       xlab="Number of heads desired", ylab="Probability", 
       titletxt="Probability of at least x heads given 10 coins thrown")


Mean and variance

Mean: $\epsilon[r]$ $$\epsilon[r] \equiv \sum_{r=0}^{N} P(r|f,N)r$$

Variance: $var[r]$ $$var[r] \equiv \epsilon [(r- \epsilon [r])^2]$$


In [142]:
def mean(N, r, f):
    a1 = [bin(N, x, f) for x in range(0, N)] # probability of each r
    a2 = sum(a1) 
    a3 = a2 * r # add prob of each r together, and multiply with r
    return(a3)

def variance(N, r, f):
    a1 = (r - mean(N, r, f)) ** 2
    a2 = mean(N, a1, f)
    return(a2)

Note that

$$\epsilon[x+y] = \epsilon[x] + \epsilon[y]$$$$var[x+y] = var[x]+var[y]$$

if x and y are independent (for example by coin tossing)

Thus mean of r is the sum of the means of those random variables, variance of r is the sum of their variances.

The mean number of heads in a single toss is $[f \cdot 1^2 + (1-f) \cdot 0^2] - f^2 = f - f^2 = f(1-f)$

Thus, mean and variance of f are: $$\epsilon[r] = N f$$ and $$var[r] = N f(1-f)$$


In [148]:
def mean2(N, r, f):
    return(N*f)

def variance2(N, r, f):
    return(N * f * (1-f))

In [147]:
mean2(10, 5, 0.5)


Out[147]:
5.0

In [4]:
bin(N=10, r=4, f=0.5)


Out[4]:
0.205078125

In [ ]: