🌐 EN | 🇯🇵 JP

Chapter 2: Probability Distributions

Mathematically Modeling How Data Is Generated

📖 Reading Time: 20-25 minutes 📊 Difficulty: Beginner 💻 Code Examples: 11

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.

🎯 Learning Objectives

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)$$

💡 Relationship Between the Bernoulli and Binomial Distributions

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
📝 Note: Why Are P(X=2) and P(X=3) the Same Value?

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$$

💡 Relationship Between the Exponential and Gamma Distributions

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.

💡 The 68-95-99.7 Rule
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.

💡 Statement of the Central Limit Theorem

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$).

⚠️ Conditions for Applying the Central Limit Theorem

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$$

⚠️ Difference Between the MLE Variance Estimator and the Unbiased Estimator

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()
💡 How to Read a Q-Q Plot

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.

✅ Review of Learning Objectives
🔑 Key Points

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.