Introduction
In Chapter 1, we learned descriptive statistics for summarizing the characteristics of data, along with the basics of probability theory for handling uncertainty (the axioms of probability, conditional probability, Bayes' theorem, expected value, and variance). Building on that foundation, this chapter takes a closer look at the Probability Distribution, which describes the correspondence between the values a random variable can take and their probabilities.
In machine learning, the model and the evaluation method you use depend on which probability distribution you assume generated the observed data. For example, the Bernoulli distribution suits binary data such as the heads or tails of a coin, the Poisson distribution suits the number of events that occur within a given period, and the normal distribution suits many natural phenomena and measurement errors. Learning to select an appropriate probability distribution is an important first step in data analysis and model design.
- Understand the characteristics of major probability distributions
- Be able to select an appropriate probability distribution
- Understand the meaning and importance of the Central Limit Theorem
- Manipulate probability distributions using SciPy
1. Discrete Probability Distributions
A Discrete Probability Distribution is a distribution for cases where the values a random variable can take are discrete (countable values such as 0, 1, 2, ...). The function that gives the probability corresponding to each value is called the Probability Mass Function (PMF).
1.1 Bernoulli Distribution
The Bernoulli Distribution is the simplest probability distribution, taking only two outcomes: success (1) or failure (0). It is used to model binary outcomes, such as a coin flip or whether a click occurred.
If the probability of success is $p$, the PMF can be written as:
$$P(X=k) = p^k (1-p)^{1-k}, \quad k \in \{0, 1\}$$
Expected value and variance:
$$E[X] = p, \qquad \text{Var}(X) = p(1-p)$$
1.2 Binomial Distribution
The Binomial Distribution is the distribution of the number of successes when a Bernoulli trial with success probability $p$ is repeated independently $n$ times. It is used in situations like "the number of heads when a coin is flipped 10 times."
$$P(X=k) = \binom{n}{k} p^k (1-p)^{n-k}, \quad k = 0, 1, \ldots, n$$
Here, $\binom{n}{k} = \frac{n!}{k!(n-k)!}$ is the binomial coefficient (the number of ways to choose $k$ items out of $n$).
Expected value and variance:
$$E[X] = np, \qquad \text{Var}(X) = np(1-p)$$
The Bernoulli distribution corresponds to the special case of the binomial distribution where the number of trials is $n=1$.
1.3 Poisson Distribution
The Poisson Distribution models the number of rare events occurring per unit of time or unit of area. It is widely used to count "how many times something happens within a given period," such as the number of calls arriving at a call center, the number of visits to a website, or the number of defective items produced at a factory.
Letting the average rate of occurrence be $\lambda$ (lambda):
$$P(X=k) = \frac{\lambda^k e^{-\lambda}}{k!}, \quad k = 0, 1, 2, \ldots$$
Expected value and variance:
$$E[X] = \lambda, \qquad \text{Var}(X) = \lambda$$
A characteristic property of the Poisson distribution is that its expected value and variance are equal, both taking the value $\lambda$. It is also known that if, in the binomial distribution, the number of trials $n$ is made large and the success probability $p$ is made small while keeping $np=\lambda$ constant, the binomial distribution approaches the Poisson distribution.
1.4 Python Implementation
Using SciPy's scipy.stats module, you can easily compute the PMF, expected value, and variance of each distribution.
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
# --- Bernoulli distribution Bernoulli(p) ---
p = 0.3
bernoulli = stats.bernoulli(p)
print("=== Bernoulli Distribution Bernoulli(p=0.3) ===")
print(f"P(X=0) = {bernoulli.pmf(0):.4f}")
print(f"P(X=1) = {bernoulli.pmf(1):.4f}")
print(f"Expected value E[X] = {bernoulli.mean():.4f}")
print(f"Variance Var(X) = {bernoulli.var():.4f}")
# --- Binomial distribution Binomial(n, p) ---
n, p = 10, 0.3
binomial = stats.binom(n, p)
print("\n=== Binomial Distribution Binomial(n=10, p=0.3) ===")
for k in [0, 3, 5, 10]:
print(f"P(X={k}) = {binomial.pmf(k):.4f}")
print(f"Expected value E[X] = {binomial.mean():.4f}")
print(f"Variance Var(X) = {binomial.var():.4f}")
# --- Poisson distribution Poisson(lambda) ---
lam = 3
poisson = stats.poisson(lam)
print("\n=== Poisson Distribution Poisson(λ=3) ===")
for k in [0, 2, 3, 6]:
print(f"P(X={k}) = {poisson.pmf(k):.4f}")
print(f"Expected value E[X] = {poisson.mean():.4f}")
print(f"Variance Var(X) = {poisson.var():.4f}")
# --- Visualization ---
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
k_bern = [0, 1]
axes[0].bar(k_bern, bernoulli.pmf(k_bern), color='skyblue', edgecolor='black')
axes[0].set_title('Bernoulli Distribution (p=0.3)')
axes[0].set_xlabel('k')
axes[0].set_ylabel('Probability')
axes[0].set_xticks(k_bern)
k_binom = np.arange(0, n + 1)
axes[1].bar(k_binom, binomial.pmf(k_binom), color='lightgreen', edgecolor='black')
axes[1].set_title('Binomial Distribution (n=10, p=0.3)')
axes[1].set_xlabel('k')
axes[1].set_ylabel('Probability')
k_poisson = np.arange(0, 12)
axes[2].bar(k_poisson, poisson.pmf(k_poisson), color='salmon', edgecolor='black')
axes[2].set_title('Poisson Distribution (λ=3)')
axes[2].set_xlabel('k')
axes[2].set_ylabel('Probability')
plt.tight_layout()
plt.show()
Execution Result:
=== Bernoulli Distribution Bernoulli(p=0.3) ===
P(X=0) = 0.7000
P(X=1) = 0.3000
Expected value E[X] = 0.3000
Variance Var(X) = 0.2100
=== Binomial Distribution Binomial(n=10, p=0.3) ===
P(X=0) = 0.0282
P(X=3) = 0.2668
P(X=5) = 0.1029
P(X=10) = 0.0000
Expected value E[X] = 3.0000
Variance Var(X) = 2.1000
=== Poisson Distribution Poisson(λ=3) ===
P(X=0) = 0.0498
P(X=2) = 0.2240
P(X=3) = 0.2240
P(X=6) = 0.0504
Expected value E[X] = 3.0000
Variance Var(X) = 3.0000
In the Poisson distribution, when $\lambda$ is an integer, the property $P(X=\lambda-1) = P(X=\lambda)$ holds. Since $\lambda=3$ here, $P(X=2)$ and $P(X=3)$ turn out to be equal. This reflects the fact that the mode of the Poisson distribution lies near $\lambda$.
2. Continuous Probability Distributions
A Continuous Probability Distribution is a distribution for cases where the random variable takes continuous (real-number) values. Since the probability of any single exact value is always 0, we instead use the Probability Density Function (PDF) $f(x)$ to express, via integration, the probability that a value falls within a given interval.
$$P(a \leq X \leq b) = \int_{a}^{b} f(x)\, dx$$
2.1 Normal Distribution
The Normal Distribution, also called the Gaussian Distribution, is a bell-shaped distribution that appears in many natural phenomena, such as height, measurement error, and test scores. It is characterized by two parameters: the mean $\mu$ and the variance $\sigma^2$.
$$f(x) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$$
Expected value and variance:
$$E[X] = \mu, \qquad \text{Var}(X) = \sigma^2$$
2.2 Exponential Distribution
The Exponential Distribution models the "waiting time" from the occurrence of one event until the next occurs. It is used for things like the interval between machine failures or the time until a customer arrives. There is a relationship in which the waiting time between events described by a Poisson distribution follows an exponential distribution.
For a rate $\lambda \gt 0$:
$$f(x) = \lambda e^{-\lambda x}, \quad x \geq 0$$
Expected value and variance:
$$E[X] = \frac{1}{\lambda}, \qquad \text{Var}(X) = \frac{1}{\lambda^2}$$
The exponential distribution has the characteristic property of memorylessness (Memorylessness). The condition "having already waited $s$ time units" does not affect the probability of waiting an additional $t$ time units:
$$P(X \gt s+t \mid X \gt s) = P(X \gt t)$$
2.3 Gamma Distribution
The Gamma Distribution generalizes the exponential distribution, representing "the waiting time until $k$ independent, exponentially distributed events have occurred." It has a shape parameter $k$ (shape) and a scale parameter $\theta$ (scale).
$$f(x) = \frac{1}{\Gamma(k)\theta^{k}} x^{k-1} e^{-x/\theta}, \quad x \geq 0$$
Here, $\Gamma(k)$ is the gamma function. Expected value and variance:
$$E[X] = k\theta, \qquad \text{Var}(X) = k\theta^2$$
Setting the shape parameter $k=1$ in the gamma distribution recovers the exponential distribution. The gamma distribution can be understood as a generalization of the exponential distribution to "multiple occurrences."
2.4 Python Implementation
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
# --- Normal distribution Normal(mu, sigma) ---
mu, sigma = 0, 1
normal = stats.norm(mu, sigma)
# --- Exponential distribution Exponential(lambda) ---
rate = 1.0 # Rate λ (in scipy, specify scale = 1/λ)
exponential = stats.expon(scale=1 / rate)
# --- Gamma distribution Gamma(k, theta) ---
shape, scale = 2.0, 2.0 # Shape parameter k, scale parameter θ
gamma_dist = stats.gamma(a=shape, scale=scale)
print("=== PDF Values and Summary Statistics for Each Distribution ===")
print(f"Normal distribution N(0,1): f(0) = {normal.pdf(0):.4f}, "
f"Mean = {normal.mean():.4f}, Variance = {normal.var():.4f}")
print(f"Exponential distribution Exp(λ=1): f(1) = {exponential.pdf(1):.4f}, "
f"Mean = {exponential.mean():.4f}, Variance = {exponential.var():.4f}")
print(f"Gamma distribution Gamma(k=2,θ=2): f(2) = {gamma_dist.pdf(2):.4f}, "
f"Mean = {gamma_dist.mean():.4f}, Variance = {gamma_dist.var():.4f}")
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
x_norm = np.linspace(-4, 4, 300)
axes[0].plot(x_norm, normal.pdf(x_norm), color='steelblue', linewidth=2)
axes[0].fill_between(x_norm, normal.pdf(x_norm), alpha=0.2, color='steelblue')
axes[0].set_title('Normal Distribution N(0, 1)')
axes[0].set_xlabel('x')
axes[0].set_ylabel('Probability Density f(x)')
x_exp = np.linspace(0, 6, 300)
axes[1].plot(x_exp, exponential.pdf(x_exp), color='darkorange', linewidth=2)
axes[1].fill_between(x_exp, exponential.pdf(x_exp), alpha=0.2, color='darkorange')
axes[1].set_title('Exponential Distribution Exp(λ=1)')
axes[1].set_xlabel('x')
axes[1].set_ylabel('Probability Density f(x)')
x_gamma = np.linspace(0, 16, 300)
axes[2].plot(x_gamma, gamma_dist.pdf(x_gamma), color='seagreen', linewidth=2)
axes[2].fill_between(x_gamma, gamma_dist.pdf(x_gamma), alpha=0.2, color='seagreen')
axes[2].set_title('Gamma Distribution Gamma(k=2, θ=2)')
axes[2].set_xlabel('x')
axes[2].set_ylabel('Probability Density f(x)')
plt.tight_layout()
plt.show()
Execution Result:
=== PDF Values and Summary Statistics for Each Distribution ===
Normal distribution N(0,1): f(0) = 0.3989, Mean = 0.0000, Variance = 1.0000
Exponential distribution Exp(λ=1): f(1) = 0.3679, Mean = 1.0000, Variance = 1.0000
Gamma distribution Gamma(k=2,θ=2): f(2) = 0.1839, Mean = 4.0000, Variance = 8.0000
3. Properties of the Normal Distribution and the Central Limit Theorem
3.1 The Standard Normal Distribution and the 68-95-99.7 Rule
A normal distribution with mean $\mu=0$ and variance $\sigma^2=1$ is called the Standard Normal Distribution. Any random variable $X$ following a normal distribution can be converted to the standard normal distribution through the following standardization.
$$z = \frac{x - \mu}{\sigma}$$
This value of $z$ is called the z-score, and it represents "how many standard deviations away from the mean" a value is.
The normal distribution has an empirical rule describing how much of the data falls within a given number of standard deviations of the mean.
- About 68.3% of the data falls within the range $\mu \pm 1\sigma$
- About 95.4% of the data falls within the range $\mu \pm 2\sigma$
- About 99.7% of the data falls within the range $\mu \pm 3\sigma$
import numpy as np
from scipy import stats
mu, sigma = 170, 8 # Example: adult male height (mean 170cm, standard deviation 8cm)
normal = stats.norm(mu, sigma)
print("=== Verifying the 68-95-99.7 Rule ===")
for k in [1, 2, 3]:
lower, upper = mu - k * sigma, mu + k * sigma
prob = normal.cdf(upper) - normal.cdf(lower)
print(f"Probability of falling within μ ± {k}σ range [{lower}, {upper}]: {prob:.4f}")
# Compute standardization (z-score)
x = 186 # A person's height
z = (x - mu) / sigma
print(f"\nStandardized score for height {x}cm: z = {z:.4f}")
print(f"Cumulative probability under the standard normal distribution (probability of being at or below this height): {stats.norm.cdf(z):.4f}")
Execution Result:
=== Verifying the 68-95-99.7 Rule ===
Probability of falling within μ ± 1σ range [162, 178]: 0.6827
Probability of falling within μ ± 2σ range [154, 186]: 0.9545
Probability of falling within μ ± 3σ range [146, 194]: 0.9973
Standardized score for height 186cm: z = 2.0000
Cumulative probability under the standard normal distribution (probability of being at or below this height): 0.9772
3.2 The Central Limit Theorem
The Central Limit Theorem (CLT) is one of the most important theorems in statistics. It states the following.
Regardless of the shape of the original distribution, the sample mean $\bar{X}_n$ of independent and identically distributed random variables $X_1, X_2, \ldots, X_n$ (with mean $\mu$ and variance $\sigma^2$) approaches a normal distribution as the sample size $n$ grows large.
$$\bar{X}_n = \frac{1}{n}\sum_{i=1}^{n} X_i \quad \xrightarrow{\ n \to \infty\ } \quad N\left(\mu, \frac{\sigma^2}{n}\right)$$
In standardized form, this becomes:
$$\frac{\bar{X}_n - \mu}{\sigma/\sqrt{n}} \xrightarrow{d} N(0, 1)$$
The reason the Central Limit Theorem is so important is that, regardless of whether the underlying population distribution is uniform (like the roll of a die) or skewed (like a binomial distribution), the sample mean approaches a normal distribution once a sufficient number of samples have been collected. This makes it possible to perform calculations that assume normality in many statistical inference procedures, such as constructing confidence intervals or conducting hypothesis tests.
3.3 Simulation in Python
Let's verify the Central Limit Theorem with a dice-rolling experiment. A single die roll follows a uniform distribution from 1 to 6, but the average of multiple rolls should approach a normal distribution.
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
def simulate_sample_means(sample_size, n_experiments=5000):
"""Repeat an experiment of rolling sample_size dice n_experiments times,
and return the sample mean for each experiment"""
rolls = np.random.randint(1, 7, size=(n_experiments, sample_size))
return rolls.mean(axis=1)
# Theoretical mean and variance of a single die
dice_mean = np.mean(np.arange(1, 7)) # 3.5
dice_var = np.var(np.arange(1, 7)) # 35/12 ≈ 2.9167
sample_sizes = [1, 2, 5, 30]
fig, axes = plt.subplots(1, 4, figsize=(18, 4))
for ax, n in zip(axes, sample_sizes):
means = simulate_sample_means(n)
ax.hist(means, bins=30, density=True, color='mediumpurple',
edgecolor='black', alpha=0.7)
ax.set_title(f'n = {n}')
ax.set_xlabel('Sample Mean')
ax.set_ylabel('Density')
print(f"n={n}: Theoretical variance (σ²/n)={dice_var/n:.4f}, "
f"Simulated variance={means.var():.4f}")
plt.suptitle('Central Limit Theorem: Distribution of Dice Sample Means', fontsize=14)
plt.tight_layout()
plt.show()
In this simulation, when the sample size is $n=1$, the distribution is simply the uniform distribution of the die outcomes themselves (each value from 1 to 6 with roughly equal probability), but as $n$ increases, the histogram gradually approaches a bell shape, and by around $n=30$ it becomes almost indistinguishable from a normal distribution. You can also confirm that the variance of the sample mean approaches the theoretical value $\sigma^2/n$ (the variance of a single die, $35/12$, divided by $n$).
The Central Limit Theorem provides a better approximation when the sample size $n$ is sufficiently large. If the original distribution is already close to normal, the approximation is good even for small $n$, but for extremely skewed distributions, a larger $n$ is required. As a rule of thumb, $n \geq 30$ is often used, but this is not an absolute standard.
4. Estimating Parameters of Probability Distributions
4.1 The Idea Behind Maximum Likelihood Estimation
In real data analysis, even when we know which probability distribution generated the observed data, the parameters of that distribution (such as $\mu, \sigma$ for the normal distribution, or $\lambda$ for the Poisson distribution) are usually unknown. Maximum Likelihood Estimation (MLE) is a representative technique for finding the most "plausible" parameters from observed data.
Assuming that the observed data $x_1, x_2, \ldots, x_n$ were drawn independently from the same distribution, the joint probability of observing this data given parameter $\theta$ (the Likelihood Function) can be written as:
$$L(\theta) = \prod_{i=1}^{n} f(x_i; \theta)$$
To simplify the calculation, it is common to maximize the Log-Likelihood, the logarithm of the likelihood, instead.
$$\ell(\theta) = \sum_{i=1}^{n} \ln f(x_i; \theta)$$
The maximum likelihood estimator $\hat{\theta}$ is defined as the parameter that maximizes this log-likelihood.
$$\hat{\theta} = \underset{\theta}{\operatorname{argmax}}\ \ell(\theta)$$
4.2 Estimating Parameters of the Normal Distribution
For the normal distribution, maximizing the log-likelihood can be shown mathematically to yield the following, intuitively natural, estimators.
$$\hat{\mu} = \bar{x} = \frac{1}{n}\sum_{i=1}^{n} x_i, \qquad \hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^{n} (x_i - \bar{x})^2$$
The unbiased estimator of sample variance we learned in Chapter 1 divides by $n-1$, but the maximum likelihood estimator of variance divides by $n$. As a result, the MLE variance estimator is known to be slightly biased downward (an underestimate). As the sample size $n$ grows, this difference becomes negligibly small.
In SciPy, the scipy.stats.<distribution name>.fit() method lets you automatically perform maximum likelihood estimation from data.
import numpy as np
from scipy import stats
np.random.seed(0)
# True parameters (normally unknown, but set here deliberately for verification)
true_mu, true_sigma = 50, 5
data = np.random.normal(true_mu, true_sigma, size=200)
# Manually compute the maximum likelihood estimators
mle_mu = np.mean(data)
mle_sigma2 = np.mean((data - mle_mu) ** 2) # MLE divides by n
mle_sigma = np.sqrt(mle_sigma2)
# Estimation via scipy.stats.norm.fit() (performs maximum likelihood estimation internally)
fit_mu, fit_sigma = stats.norm.fit(data)
print("=== Estimating Normal Distribution Parameters (Maximum Likelihood) ===")
print(f"True parameters: μ = {true_mu}, σ = {true_sigma}")
print(f"Manual MLE estimates: μ_hat = {mle_mu:.4f}, σ_hat = {mle_sigma:.4f}")
print(f"scipy.stats.norm.fit: μ_hat = {fit_mu:.4f}, σ_hat = {fit_sigma:.4f}")
The manually computed estimates and the results from scipy.stats.norm.fit() agree (depending on the random seed, some deviation from the true parameters will occur due to the finite sample size). As the sample size $n$ increases, the estimates approach the true parameters.
4.3 Estimating Parameters of the Binomial Distribution
For the success probability $p$ of the Bernoulli and binomial distributions as well, the maximum likelihood estimator can be shown to equal the observed proportion of successes.
$$\hat{p} = \frac{\text{Number of successes}}{\text{Number of trials}}$$
import numpy as np
from scipy import stats
# Experiment: flipping 200 coins, 140 came up heads
n_trials = 200
n_success = 140
# Maximum likelihood estimate: p_hat = number of successes / number of trials
p_mle = n_success / n_trials
print("=== Estimating Binomial Distribution Parameters (Maximum Likelihood) ===")
print(f"Number of trials: {n_trials}, Number of successes: {n_success}")
print(f"Maximum likelihood estimate p_hat = {p_mle:.4f}")
# Build a binomial distribution using the estimated p and check its summary statistics
binom_est = stats.binom(n_trials, p_mle)
print(f"Expected value under the estimated distribution: {binom_est.mean():.4f}")
print(f"Standard deviation under the estimated distribution: {binom_est.std():.4f}")
Execution Result:
=== Estimating Binomial Distribution Parameters (Maximum Likelihood) ===
Number of trials: 200, Number of successes: 140
Maximum likelihood estimate p_hat = 0.7000
Expected value under the estimated distribution: 140.0000
Standard deviation under the estimated distribution: 6.4807
5. Visualizing and Simulating Probability Distributions
5.1 Overlaying Histograms with Theoretical Distributions
A basic way to check whether observed data follows a particular probability distribution is to overlay the PDF of a theoretical distribution, with parameters estimated from the data, onto a histogram of that data.
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(1)
# Example: reaction time data (assumed to follow a gamma distribution)
data = np.random.gamma(shape=3.0, scale=1.5, size=500)
# Estimate the gamma distribution parameters from the data (maximum likelihood)
# floc=0 fixes the location parameter to 0, estimating only shape and scale
fit_shape, fit_loc, fit_scale = stats.gamma.fit(data, floc=0)
print("=== Estimating Gamma Distribution Parameters ===")
print(f"Estimated shape parameter k = {fit_shape:.4f}")
print(f"Estimated scale parameter θ = {fit_scale:.4f}")
# Overlay the histogram with the estimated distribution
x = np.linspace(0, data.max(), 300)
fitted_pdf = stats.gamma.pdf(x, a=fit_shape, scale=fit_scale)
plt.figure(figsize=(10, 6))
plt.hist(data, bins=30, density=True, alpha=0.6, color='lightsteelblue',
edgecolor='black', label='Observed data')
plt.plot(x, fitted_pdf, color='crimson', linewidth=2, label='Fitted gamma distribution')
plt.xlabel('Value')
plt.ylabel('Density')
plt.title('Histogram Overlaid with the Fitted Distribution')
plt.legend()
plt.grid(alpha=0.3)
plt.show()
The estimated shape parameter $k$ and scale parameter $\theta$ come out close to the true values used to generate the data ($k=3.0$, $\theta=1.5$), and you can confirm that the fitted curve closely matches the shape of the histogram.
5.2 Checking Distribution Fit with a Q-Q Plot
A Q-Q Plot (Quantile-Quantile Plot) is a technique for comparing the quantiles of observed data against the quantiles of an assumed theoretical distribution as a scatter plot. If the data follows the theoretical distribution, the points will fall approximately on a straight line.
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(2)
# Compare normal-distribution data with non-normal (exponential-distribution) data
normal_data = np.random.normal(0, 1, 300)
exp_data = np.random.exponential(scale=1.0, size=300)
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
stats.probplot(normal_data, dist="norm", plot=axes[0])
axes[0].set_title('Q-Q Plot: Normal Distribution Data')
stats.probplot(exp_data, dist="norm", plot=axes[1])
axes[1].set_title('Q-Q Plot: Exponential Distribution Data (Compared to Normal)')
plt.tight_layout()
plt.show()
- Points lying roughly along the straight line: the data closely follows the assumed distribution (here, the normal distribution)
- Points systematically deviating from the line (especially curving at both ends): the data likely does not follow the assumed distribution
- The Q-Q plot for normal-distribution data lies close to a straight line, while the Q-Q plot for the right-skewed exponential-distribution data deviates substantially from the line at both ends, especially in the upper right
6. Summary and Next Steps
In this chapter, we learned about probability distributions that frequently appear in machine learning, along with methods for characterizing them from data.
- Understand the characteristics of major probability distributions: We learned the PMF/PDF, expected value, and variance of the Bernoulli, binomial, and Poisson distributions (discrete) and the normal, exponential, and gamma distributions (continuous)
- Be able to select an appropriate probability distribution: We confirmed that the Bernoulli distribution suits binary data, the Poisson and binomial distributions suit count data, the normal distribution suits continuous measurements, and the gamma and exponential distributions suit waiting times
- Understand the meaning and importance of the Central Limit Theorem: We used a dice-rolling simulation to confirm that, regardless of the shape of the original distribution, the sample mean approaches a normal distribution
- Manipulate probability distributions using SciPy: We implemented PMF/PDF calculations, parameter estimation (
fit()), visualization, and fit-checking with Q-Q plots usingscipy.stats
- Discrete distributions express probability with a PMF, and continuous distributions with a PDF
- The Poisson distribution has equal expected value and variance, and the exponential distribution is memoryless
- By the Central Limit Theorem, the sample mean approaches a normal distribution regardless of the shape of the original distribution
- Maximum likelihood estimation finds the parameters that make the observed data most probable
- Overlaying histograms with fitted curves and using Q-Q plots let you visually check whether an assumed distribution fits the data
Next Steps
In the next chapter, building on the knowledge of probability distributions from this chapter, we will learn about statistical estimation and hypothesis testing, which infer properties of a population from a sample. We will cover techniques essential for data-driven decision making, such as constructing confidence intervals and interpreting p-values.
Practice Problems
Problem 1: Application of the Poisson Distribution
A call center receives, on average, 4 calls per hour. Using the Poisson distribution, write code to compute the probability that at most 2 calls are received in one hour.
from scipy import stats
lam = 4
poisson = stats.poisson(lam)
# Sum of the probabilities of 2 or fewer calls (0, 1, or 2). Using the cumulative distribution function (CDF) is efficient
prob = poisson.cdf(2)
print(f"Probability that at most 2 calls are received in one hour: {prob:.4f}")
Answer: The probability of receiving at most 2 calls in one hour is about 0.2381 (23.81%). By using the Poisson distribution's cumulative distribution function cdf(), you can find $P(X=0)+P(X=1)+P(X=2)$ without computing each term individually.
Problem 2: Simulating the Central Limit Theorem
Assume the population follows a Uniform(0, 10) distribution. Write code that draws 1000 samples of size 50, visualizes the distribution of the sample means as a histogram, and compares the theoretical mean and variance (population mean, population variance / n) with the simulation results.
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
population_mean = 5.0 # Theoretical mean of Uniform(0,10), (a+b)/2
population_var = 100 / 12 # Theoretical variance of Uniform(0,10), (b-a)^2/12
n_experiments = 1000
sample_size = 50
sample_means = np.array([
np.random.uniform(0, 10, sample_size).mean()
for _ in range(n_experiments)
])
theoretical_mean = population_mean
theoretical_var = population_var / sample_size
print(f"Theoretical mean of the sample mean: {theoretical_mean:.4f}")
print(f"Theoretical variance of the sample mean: {theoretical_var:.4f}")
print(f"Simulated mean of the sample mean: {sample_means.mean():.4f}")
print(f"Simulated variance of the sample mean: {sample_means.var():.4f}")
plt.hist(sample_means, bins=30, density=True, color='teal', alpha=0.7, edgecolor='black')
plt.xlabel('Sample Mean')
plt.ylabel('Density')
plt.title('Distribution of Sample Means from a Uniform Distribution (Central Limit Theorem)')
plt.show()
Answer: The theoretical mean of the sample mean is 5.0000, and the theoretical variance is $100/12/50 \approx 0.1667$. The simulation results come out close to these theoretical values, and with a sample size of 50, the histogram of the sample mean is nearly normal in shape even though the underlying distribution is uniform.
Problem 3: Practicing Maximum Likelihood Estimation
Assume that the observed data for the number of defective items produced per day at a factory follows a Poisson distribution. Write code to estimate the parameter $\lambda$ from the following data using maximum likelihood estimation.
Data: [2, 4, 3, 5, 2, 6, 3, 4, 1, 5]
import numpy as np
data = np.array([2, 4, 3, 5, 2, 6, 3, 4, 1, 5])
# The maximum likelihood estimator of the Poisson distribution parameter λ equals the sample mean
lambda_mle = np.mean(data)
print(f"Maximum likelihood estimate λ_hat = {lambda_mle:.4f}")
Answer: Since the sum of the data is 35 and there are 10 data points, the maximum likelihood estimate is $\hat{\lambda} = 35/10 = 3.5000$. Differentiating the Poisson distribution's log-likelihood with respect to $\lambda$ and setting it to zero mathematically shows that the maximum likelihood estimator equals the sample mean.