From Chapter 1 of Information Theory, Inference and Learning Algorithms
A coin has probability f
of coming up heads, tossed N
times. Probability of getting r
heads? Mean and variance of 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)
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")
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")
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")
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]:
In [4]:
bin(N=10, r=4, f=0.5)
Out[4]:
In [ ]: