Understanding Binomial Distribution, Permutations, Combinations, and Their Python Implementations

Binomial distribution is a basic and practical probability model that describes how many times a particular outcome occurs in a fixed number of independent trials.

Bernoulli distribution, the foundation of the binomial model, has two possible outcomes (0 and 1) with probabilities that sum to 1; its probability mass function, expectation, and variance are shown in the accompanying images.

Permutations count ordered selections of r objects from n objects, while combinations count unordered selections; the formulas for both are illustrated with images.

The following Python function PC(n, k) computes the permutation and combination numbers for given n and k :

<code>from functools import reduce

def PC(n, k):
    """Calculate and return permutation and combination numbers"""
    # Invalid input returns None
    if n <= 0 or k < 0 or n < k:
        print('Wrong Input!')
        return None
    # When k is 0, both values are 1
    if k == 0:
        return 1, 1
    series_asc = list(range(1, n+1))
    series_desc = sorted(series_asc, reverse=True)
    permutation = reduce(lambda x, y: x*y, series_desc[:k])
    perm2 = reduce(lambda x, y: x*y, series_asc[:k])
    combination = int(permutation / perm2)
    return permutation, combination</code>

Testing the function with several n and k values confirms its correctness.

A second function binominal_prob(n, k, p) computes the probability of obtaining exactly k successes in n independent Bernoulli trials with success probability p :

<code>def binominal_prob(n, k, p):
    """Calculate the binomial probability for given n, k, and p"""
    p_base = p ** k * (1 - p) ** (n - k)
    combination = PC(n, k)[1]
    return p_base * combination</code>

Using this function, probabilities for n = 10 trials are computed for all possible numbers of successes and visualized with matplotlib :

<code>%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
probs = [round(i/10,1) for i in range(1,10)]
n = 20
plt.figure(figsize=(16,6))
for p in probs:
    dist_probs = [binominal_prob(n, i, p) for i in range(n+1)]
    plt.plot(range(n+1), dist_probs, label='p={0}'.format(p))
plt.legend()
plt.title('Binominal Distributions of Different P value')
plt.savefig('binominal.jpg')
plt.show()</code>

The same calculation is reproduced with the interactive plotly library:

<code>import plotly.graph_objects as go
probs = [round(i/10,1) for i in range(1,10)]
n = 20
fig = go.Figure()
for p in probs:
    dist_probs = [binominal_prob(n, i, p) for i in range(n+1)]
    fig.add_trace(go.Scatter(x=list(range(n+1)), y=dist_probs, name='p={0}'.format(p)))
fig.show()</code>

Finally, a maximum‑likelihood estimation example shows how to estimate the success probability p of a biased coin from four flips that resulted in three heads; differentiating the likelihood function yields the optimal estimate p = 0.75 , which maximizes the probability of the observed outcome.