Probability and Statistics Practice

The Probability and Statistics chapter established the theoretical foundations of probability and statistics. This chapter puts that knowledge into practice, using NumPy to implement core probability and statistical computations. Through hands-on coding, you will not only deepen your understanding of the concepts but also acquire practical data analysis skills.

Distribution Sampling

The Uniform Distribution is the simplest continuous probability distribution, with a probability density function that is constant throughout its defined interval. For a uniform distribution over $[a, b]$, the probability that a variable $X$ falls within any subinterval $[c, d]$ is proportional to the subinterval length: $P(c \leq X \leq d) = \frac{d-c}{b-a}$. Uniform distributions are commonly used in practice to simulate "unbiased" random selection, such as random sample selection or random task assignment.

The Normal Distribution (also called the Gaussian distribution) is the most important distribution in statistics, with a symmetric bell-shaped probability density function. Many phenomena in nature and human society approximately follow a normal distribution, such as human height, exam scores, measurement errors, etc. The standard normal distribution $N(0, 1)$ has a probability density function: $f(x) = \frac{1}{\sqrt{2\pi}} e^{-x^2/2}$, characterized by a mean (expectation) of 0 and a standard deviation of 1. The curve peaks at the mean and gradually decreases toward both sides, with approximately 68% of the data falling within $\pm 1$ standard deviation and about 95% within $\pm 2$ standard deviations.

Sampling from a probability distribution means generating random numbers according to that distribution's probability density function — the higher the probability of a value, the more likely it is to be generated. Sampling is a fundamental operation for Monte Carlo simulation, randomized algorithms, statistical inference, and other techniques. NumPy provides concise functions for sampling from these two most basic distributions:

• np.random.rand(n): Samples from a uniform distribution over $[0, 1)$, returning $n$ samples.
• np.random.randn(n): Samples from the standard normal distribution $N(0, 1)$, returning $n$ samples.
• np.random.uniform(low, high, n): Samples from a uniform distribution over an arbitrary interval $[low, high]$.
• np.random.normal(mean, std, n): Samples from a normal distribution $N(\mu, \sigma^2)$ with arbitrary parameters.

The following code demonstrates the sampling process for uniform and normal distributions, comparing the histograms of the sampled data with the theoretical probability density function (PDF) via visualization to verify the correctness of sampling.

import numpy as np
import matplotlib.pyplot as plt

# Uniform distribution sampling
n = 10000
uniform_samples = np.random.rand(n)  # [0, 1) uniform distribution

# Normal distribution sampling
normal_samples = np.random.randn(n)  # Standard normal distribution

# Visualization
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Uniform distribution
axes[0].hist(uniform_samples, bins=50, density=True, alpha=0.7, color='steelblue', edgecolor='black')
axes[0].axhline(1, color='r', linestyle='--', linewidth=2, label='Theoretical PDF')
axes[0].set_xlabel('Value')
axes[0].set_ylabel('Probability Density')
axes[0].set_title(f'Uniform Distribution Sampling (n={n})')
axes[0].legend()
axes[0].grid(alpha=0.3)

# Normal distribution
axes[1].hist(normal_samples, bins=50, density=True, alpha=0.7, color='steelblue', edgecolor='black')

# Theoretical PDF
x = np.linspace(-4, 4, 100)
pdf = 1 / np.sqrt(2 * np.pi) * np.exp(-x**2 / 2)
axes[1].plot(x, pdf, 'r-', linewidth=2, label='Theoretical PDF')
axes[1].set_xlabel('Value')
axes[1].set_ylabel('Probability Density')
axes[1].set_title(f'Standard Normal Distribution Sampling (n={n})')
axes[1].legend()
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.show()
plt.close()

Click Run to execute code. Click the code area to edit.

Random Seeds and Reproducibility

A Random Seed is a value that controls the initial state of a random number generator. "Randomness" in computers is actually pseudo-random, generated by deterministic algorithms. Given the same seed, the algorithm will produce the exact same sequence of random numbers. In machine learning experiments, scientific computing, and data analysis, Reproducibility is a fundamental requirement — others should be able to repeat your experiments and obtain the same results. This is indispensable for result verification, error diagnosis, and model comparison.

NumPy's random number generator maintains an internal state. Setting the random seed via np.random.seed(n) to an integer $n$ ensures that the same seed produces the same sequence across different runs, making all subsequent random operations deterministic. Commonly used seed values such as 42, 0, etc., have become community conventions (they have no special meaning — they are simply easy to remember). The following code compares the difference in random number generation without and with a seed setting, intuitively demonstrating how random seeds control reproducibility.

import numpy as np

# Without a seed: different results each time
print("Without setting a random seed:")
for i in range(3):
    samples = np.random.rand(5)
    print(f"  Run {i+1}: {samples}")

print()

# With a seed: same results each time
print("With random seed (seed=42):")
for i in range(3):
    np.random.seed(42)
    samples = np.random.rand(5)
    print(f"  Run {i+1}: {samples}")

Click Run to execute code. Click the code area to edit.

Random Selection and Shuffling

Random Selection is the operation of drawing elements from a given set according to specified rules. Depending on whether repeated draws are allowed, it can be classified as sampling with replacement (the same element can be selected multiple times) or sampling without replacement (each element can be selected at most once). Weighted selection allows assigning a probability weight to each element, which in practice can be used to simulate decision scenarios where different options have different probabilities, such as diversity sampling in recommendation systems or stratified sampling in surveys.

Random Shuffling is the random rearrangement of the order of elements in a sequence, such that each element has an equal probability of appearing at any position. Shuffling is extremely common in data preprocessing: random partitioning of training and test sets, random ordering in data augmentation, random grouping in experimental design, and other scenarios require breaking the original order of data to avoid bias.

NumPy provides the following functions to implement these operations:

• np.random.choice(array, size, replace): Randomly selects elements from an array; replace controls whether duplicates are allowed
• np.random.choice(array, size, p): Weighted selection; p specifies the probability of each element being selected (must have the same length as the array and sum to 1)
• np.random.shuffle(array): Shuffles the array in-place, directly modifying the original array

The following code demonstrates three modes of random selection (with replacement, without replacement, and weighted) as well as random shuffling.

import numpy as np
# Random selection
data = np.array(['apple', 'banana', 'orange', 'grape', 'watermelon'])

# Selection with replacement
choices_with_replacement = np.random.choice(data, size=10, replace=True)
print("Selection with replacement:", choices_with_replacement)

# Selection without replacement
choices_without_replacement = np.random.choice(data, size=3, replace=False)
print("Selection without replacement:", choices_without_replacement)

# Weighted selection
weights = [0.4, 0.3, 0.15, 0.1, 0.05]  # Probability of each element being selected
weighted_choices = np.random.choice(data, size=10, p=weights)
print("Weighted selection:", weighted_choices)

# Random shuffling
arr = np.arange(10)
print(f"\nOriginal array: {arr}")
np.random.shuffle(arr)
print(f"After shuffling: {arr}")

Click Run to execute code. Click the code area to edit.

Distribution Visualization

By visualizing the shape of distributions, we can intuitively understand their characteristics: histograms of discrete distributions (Binomial Distribution, Poisson Distribution) appear as bars, with the height of each bar reflecting the probability of that value; histograms of continuous distributions (Exponential Distribution) need parameter adjustments to present smooth curves, with height reflecting probability density. NumPy provides a rich set of distribution sampling functions:

• np.random.binomial(n, p, size): Samples from a binomial distribution $B(n, p)$
• np.random.poisson(lam, size): Samples from a Poisson distribution $Poisson(\lambda)$
• np.random.exponential(scale, size): Samples from an exponential distribution $Exp(scale)$, where scale is the average waiting time

The following code demonstrates sampling and visualization of three distributions, using np.histogram to compute raw histogram data.

import numpy as np
import matplotlib.pyplot as plt

n = 10000 # Generate data from multiple distributions
# Binomial distribution
binomial = np.random.binomial(n=20, p=0.3, size=n)
# Poisson distribution
poisson = np.random.poisson(lam=5, size=n)
# Exponential distribution
exponential = np.random.exponential(scale=2, size=n)

# Visualization
fig, axes = plt.subplots(1, 3, figsize=(14, 4))

# Binomial distribution (discrete)
axes[0].hist(binomial, bins=range(0, 21), density=True, color='steelblue', edgecolor='black', alpha=0.7)
axes[0].set_xlabel('Value')
axes[0].set_ylabel('Probability')
axes[0].set_title('Binomial Distribution B(20, 0.3)')
axes[0].grid(alpha=0.3, axis='y')

# Poisson distribution (discrete)
axes[1].hist(poisson, bins=range(0, 20), density=True, color='steelblue', edgecolor='black', alpha=0.7)
axes[1].set_xlabel('Value')
axes[1].set_ylabel('Probability')
axes[1].set_title('Poisson Distribution Poisson(5)')
axes[1].grid(alpha=0.3, axis='y')

# Exponential distribution (continuous)
axes[2].hist(exponential, bins=50, density=True, color='steelblue', edgecolor='black', alpha=0.7)
axes[2].set_xlabel('Value')
axes[2].set_ylabel('Probability Density')
axes[2].set_title('Exponential Distribution (scale=2)')
axes[2].grid(alpha=0.3)
plt.tight_layout()
plt.show()
plt.close()

# Use np.histogram to compute histogram data
hist, bin_edges = np.histogram(binomial, bins=range(0, 22), density=True)
print("Binomial distribution histogram data:")
print(f"  Bin edges: {bin_edges[:6]}...")  # Show only the first few
print(f"  Probabilities: {hist[:5]}...")

Click Run to execute code. Click the code area to edit.

Computing Statistics

Descriptive Statistics are numerical summaries of basic characteristics of a dataset, divided into three categories: central tendency (mean, median), dispersion (variance, standard deviation, range), and distribution shape (quantiles, skewness, kurtosis). These statistics do not rely on distributional assumptions and provide a "first glance" understanding of the data, helping to quickly identify data characteristics, detect outliers, and compare different datasets. NumPy provides efficient functions for computing statistics:

• np.mean(data) / np.median(data): Compute mean and median
• np.var(data) / np.std(data): Compute variance and standard deviation
• np.min(data) / np.max(data) / np.ptp(data): Minimum, maximum, range
• np.percentile(data, p): Compute the $p$-th percentile
• np.cov(x, y) / np.corrcoef(x, y): Covariance and correlation coefficient

The following code uses simulated exam score data as an example, demonstrating the computation of various descriptive statistics and manually implementing the formulas for skewness and kurtosis.

import numpy as np

data = np.random.normal(100, 15, 1000)  # Simulated exam scores

print("=== Descriptive Statistics ===")
print(f"Sample size: {len(data)}")
print(f"Minimum: {np.min(data):.2f}")
print(f"Maximum: {np.max(data):.2f}")
print(f"Range: {np.ptp(data):.2f}")  # peak-to-peak
print()
print(f"Mean: {np.mean(data):.2f}")
print(f"Median: {np.median(data):.2f}")
print()
print(f"Variance: {np.var(data):.2f}")
print(f"Standard deviation: {np.std(data):.2f}")
print()

# Percentiles
percentiles = [25, 50, 75]
for p in percentiles:
    print(f"{p}th percentile: {np.percentile(data, p):.2f}")

# Interquartile range
q1, q3 = np.percentile(data, [25, 75])
iqr = q3 - q1
print(f"\nInterquartile Range (IQR): {iqr:.2f}")

# Skewness and kurtosis (manual computation)
mean = np.mean(data)
std = np.std(data)
skewness = np.mean(((data - mean) / std) ** 3)
kurtosis = np.mean(((data - mean) / std) ** 4) - 3

print(f"\nSkewness: {skewness:.4f} (normal distribution = 0)")
print(f"Kurtosis: {kurtosis:.4f} (normal distribution = 0)")

Click Run to execute code. Click the code area to edit.

Covariance and Correlation Coefficient

Covariance measures the degree to which two random variables change together. For variables $X$ and $Y$, covariance is defined as $Cov(X, Y) = E[(X - E[X])(Y - E[Y])]$. A positive covariance indicates that the two variables tend to change in the same direction (when one increases, the other also increases), a negative value indicates opposite direction changes, and zero indicates linear independence. The magnitude of covariance is affected by the units of the variables, making it difficult to directly compare the strength of relationships across different datasets.

The Correlation Coefficient is a standardized version of covariance that eliminates the influence of units. The correlation coefficient is defined as $r = \frac{Cov(X, Y)}{\sigma_X \sigma_Y}$, with a range of $[-1, 1]$. $r = 1$ indicates a perfect positive correlation, $r = -1$ indicates a perfect negative correlation, and $r = 0$ indicates no linear correlation. The correlation coefficient is one of the most commonly used measures of relationship in data analysis, widely applied in feature selection, factor analysis, regression diagnostics, and other scenarios.

Covariance matrices and correlation matrices extend this measurement to multivariate scenarios. For $n$ variables, the element at the $i$-th row and $j$-th column represents the covariance or correlation coefficient between the $i$-th and $j$-th variables. Diagonal elements are variances (for the covariance matrix) or 1 (for the correlation matrix). NumPy provides the following functions for computing these statistics:

• np.cov(x, y): Computes the covariance of two variables
• np.cov([x1, x2, ...]): Computes the covariance matrix of multiple variables
• np.corrcoef(x, y): Computes the correlation coefficient of two variables
• np.corrcoef([x1, x2, ...]): Computes the correlation matrix of multiple variables

The following code generates three sets of data: $y$ positively correlated with $x$, and $z$ independent of $x$. It demonstrates the computation of the covariance matrix and correlation matrix, and visually displays the effect of different correlation strengths through scatter plots.

import numpy as np
import matplotlib.pyplot as plt

# Generate correlated data
n = 100
x = np.random.randn(n)
y = 0.8 * x + 0.2 * np.random.randn(n)  # y positively correlated with x
z = np.random.randn(n)  # z uncorrelated with x

# Compute covariance matrix
cov_matrix = np.cov([x, y, z])
print("Covariance matrix:")
print(np.round(cov_matrix, 3))
print()

# Compute correlation matrix
corr_matrix = np.corrcoef([x, y, z])
print("Correlation matrix:")
print(np.round(corr_matrix, 3))
print()

# Visualization
fig, axes = plt.subplots(1, 3, figsize=(12, 4))
axes[0].scatter(x, y, alpha=0.6)
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].set_title(f'x vs y (r = {corr_matrix[0,1]:.2f})')
axes[0].grid(alpha=0.3)
axes[1].scatter(x, z, alpha=0.6)
axes[1].set_xlabel('x')
axes[1].set_ylabel('z')
axes[1].set_title(f'x vs z (r = {corr_matrix[0,2]:.2f})')
axes[1].grid(alpha=0.3)
axes[2].scatter(y, z, alpha=0.6)
axes[2].set_xlabel('y')
axes[2].set_ylabel('z')
axes[2].set_title(f'y vs z (r = {corr_matrix[1,2]:.2f})')
axes[2].grid(alpha=0.3)
plt.tight_layout()
plt.show()
plt.close()

Click Run to execute code. Click the code area to edit.

Monte Carlo Method

Many practical problems are difficult to solve exactly using analytical methods — for example, integral computations where the function has no closed-form expression, probability estimation where the event combinations are too complex, optimization problems where the objective function is non-differentiable, and so on. The Monte Carlo Method provides an approximate solution path for such problems.

The Monte Carlo method is named after the Monte Carlo casino in Monaco, with randomness being the core of the method. This approach transforms complex problems into statistical aggregation of a large number of simple samples through extensive random sampling. Its working principle rests on two mathematical foundations: the Law of Large Numbers guarantees that as the sample size approaches infinity, the sample mean converges to the expected value (the true value); the Central Limit Theorem guarantees that the error distribution of convergence approximates a normal distribution, making the estimation error quantifiable and predictable. Specifically, using sample means to estimate integral values, success frequencies to estimate probabilities, and random experiment statistics to approximate analytical results — the larger the sample size, the higher the estimation accuracy.

• The following code demonstrates the Monte Carlo method for computing the integral $\int_a^b f(x) dx \approx \frac{b-a}{N} \sum_{i=1}^N f(x_i)$, where $x_i$ are uniform random samples over $[a, b]$:

  import numpy as np
  import matplotlib.pyplot as plt
  
  # Compute ∫_0^1 sin(x) dx
  # True value: 1 - cos(1) ≈ 0.4597
  def f(x):
      return np.sin(x)
  a, b = 0, 1
  true_value = 1 - np.cos(1)  # Analytical solution
  
  # Estimates with different sample sizes
  sample_sizes = [100, 1000, 10000, 100000]
  estimates = []
  
  for n in sample_sizes:
      x_samples = np.random.uniform(a, b, n)
      estimate = (b - a) * np.mean(f(x_samples))
      estimates.append(estimate)
      error = abs(estimate - true_value)
      print(f"n = {n:6d}: estimate = {estimate:.6f}, error = {error:.6f}")
  print(f"\nTrue value: {true_value:.6f}")
  
  # Visualize convergence
  n_range = np.arange(100, 10001, 100)
  convergence = []
  for n in n_range:
      x_samples = np.random.uniform(a, b, n)
      estimate = (b - a) * np.mean(f(x_samples))
      convergence.append(estimate)
  
  plt.figure(figsize=(10, 5))
  plt.plot(n_range, convergence, 'b-', alpha=0.5, label='Monte Carlo Estimate')
  plt.axhline(true_value, color='r', linestyle='--', linewidth=2, label=f'True Value = {true_value:.4f}')
  plt.xlabel('Sample Size n')
  plt.ylabel('Integral Estimate')
  plt.title('Convergence of Monte Carlo Integral Estimation')
  plt.legend()
  plt.grid(alpha=0.3)
  plt.tight_layout()
  plt.show()
  plt.close()
  

  Click Run to execute code. Click the code area to edit.

• Another classic example of the Monte Carlo method is estimating pi by randomly throwing points. Within a square of side length 2, points are thrown randomly, and the proportion falling inside the unit circle is counted. Since the ratio of the circle area $\pi$ to the square area $4$ is $\pi/4$, multiplying this proportion by 4 gives the estimate of pi.

  import numpy as np
  import matplotlib.pyplot as plt
  
  # Estimate pi using the Monte Carlo method
  # Randomly throw points in a unit square; proportion inside unit circle = pi/4
  n = 10000
  x = np.random.uniform(-1, 1, n)
  y = np.random.uniform(-1, 1, n)
  # Check if points are inside the unit circle
  inside = x**2 + y**2 <= 1
  # Estimate pi
  pi_estimate = 4 * np.sum(inside) / n
  error = abs(pi_estimate - np.pi)
  
  print(f"Number of sample points: {n}")
  print(f"Points inside the circle: {np.sum(inside)}")
  print(f"Pi estimate: {pi_estimate:.6f}")
  print(f"True value: {np.pi:.6f}")
  print(f"Error: {error:.6f}")
  
  # Visualization
  plt.figure(figsize=(8, 8))
  # Plot points
  plt.scatter(x[inside], y[inside], c='blue', s=1, alpha=0.5, label='Inside circle')
  plt.scatter(x[~inside], y[~inside], c='red', s=1, alpha=0.5, label='Outside circle')
  # Plot circle
  theta = np.linspace(0, 2*np.pi, 100)
  plt.plot(np.cos(theta), np.sin(theta), 'k-', linewidth=2)
  # Plot square
  plt.plot([-1, 1, 1, -1, -1], [-1, -1, 1, 1, -1], 'k-', linewidth=2)
  plt.xlabel('x')
  plt.ylabel('y')
  plt.title(f'Monte Carlo Estimate of pi = {pi_estimate:.4f}')
  plt.axis('equal')
  plt.legend()
  plt.grid(alpha=0.3)
  plt.tight_layout()
  plt.show()
  plt.close()
  
  # Convergence process
  sample_sizes = [100, 1000, 5000, 10000, 50000, 100000]
  estimates = []
  for n_samples in sample_sizes:
      x = np.random.uniform(-1, 1, n_samples)
      y = np.random.uniform(-1, 1, n_samples)
      inside = x**2 + y**2 <= 1
      estimates.append(4 * np.sum(inside) / n_samples)
  
  print("\nConvergence:")
  for n_samples, est in zip(sample_sizes, estimates):
      print(f"  n = {n_samples:6d}: pi = {est:.6f}, error = {abs(est - np.pi):.6f}")
  

  Click Run to execute code. Click the code area to edit.

• The Monte Carlo method can also estimate the probability of complex events. By simulating the event process through a large number of random trials, the frequency of events is used as an estimate of the probability. When the number of trials is sufficiently large, the frequency converges to the probability. The following code estimates the probability that the sum of three standard normal variables exceeds 3:

  import numpy as np
  import matplotlib.pyplot as plt
  
  # Estimate the probability that the sum of three standard normal variables exceeds 3
  def estimate_probability(n_samples=100000):
      """Estimate P(X + Y + Z > 3), where X, Y, Z ~ N(0,1)"""
      X = np.random.randn(n_samples)
      Y = np.random.randn(n_samples)
      Z = np.random.randn(n_samples)
      S = X + Y + Z
      count = np.sum(S > 3)   
      return count / n_samples
  
  # Monte Carlo estimate
  n_samples = 100000
  prob_estimate = estimate_probability(n_samples)
  
  # Theoretical value (sum S of three independent N(0,1) variables follows N(0, 3); standardized S/sqrt(3) ~ N(0,1))
  from math import erf, sqrt
  z = 3 / sqrt(3)
  prob_theory = 0.5 * (1 - erf(z / sqrt(2)))
  print(f"Estimation of P(X + Y + Z > 3)")
  print(f"  Monte Carlo estimate: {prob_estimate:.6f}")
  print(f"  Theoretical value: {prob_theory:.6f}")
  print(f"  Error: {abs(prob_estimate - prob_theory):.6f}")
  
  # Visualize the distribution of S
  n = 100000
  X = np.random.randn(n)
  Y = np.random.randn(n)
  Z = np.random.randn(n)
  S = X + Y + Z
  
  plt.figure(figsize=(10, 5))
  plt.hist(S, bins=100, density=True, alpha=0.7, color='steelblue', edgecolor='black')
  
  # Theoretical PDF
  x = np.linspace(-6, 6, 100)
  pdf = 1 / sqrt(2 * np.pi * 3) * np.exp(-x**2 / (2 * 3))
  plt.plot(x, pdf, 'r-', linewidth=2, label='Theoretical PDF: N(0, 3)')
  plt.axvline(3, color='green', linestyle='--', linewidth=2, label='Threshold = 3')
  plt.xlabel('S = X + Y + Z')
  plt.ylabel('Probability Density')
  plt.title('Distribution of the Sum of Three Standard Normal Variables')
  plt.legend()
  plt.grid(alpha=0.3)
  plt.tight_layout()
  plt.show()
  plt.close()
  

  Click Run to execute code. Click the code area to edit.

Bayesian Posterior Sampling

Bayesian Inference is a statistical inference method that updates our knowledge of unknown parameters based on observed data. Bayes' theorem combines the prior distribution of a parameter with the likelihood function of the data to obtain the posterior distribution: $P(\theta | data) = \frac{P(data | \theta) P(\theta)}{P(data)}$. Here, $P(\theta)$ is called the Prior Distribution, representing our knowledge of the parameter before observing the data; $P(data | \theta)$ is called the Likelihood Function, describing the probability of observing the data given a particular parameter value; $P(data)$ is a normalization constant ensuring the posterior distribution integrates to 1. The posterior distribution fully characterizes the uncertainty of the parameter, containing all possible values of the parameter and their relative probabilities.

Posterior Sampling is the process of drawing samples from the posterior distribution. When the posterior distribution has no analytical form or is computationally difficult, sampling methods provide a practical approach for obtaining posterior information. Through a large number of posterior samples, we can estimate the posterior mean, variance, credible intervals, and other statistics of the parameter, achieving a complete description of parameter uncertainty. Bayesian posterior sampling is widely used in machine learning for parameter estimation, model selection, predictive inference, and other scenarios.

NumPy provides the following methods to implement Bayesian posterior sampling:

• np.random.gamma(shape, scale, size): Samples from a Gamma distribution, used for Beta-Gamma relationship sampling
• np.random.binomial(n, p, size): Generates binomial observation data
• np.percentile(samples, p): Computes credible intervals from posterior samples

The following code demonstrates the complete workflow of Bayesian posterior sampling using the estimation of a coin's head probability as an example: generating observation data, computing posterior parameters, sampling from the Beta posterior distribution, and summarizing posterior information.

import numpy as np
import matplotlib.pyplot as plt

# Bayesian posterior sampling: estimating the probability of heads for a coin
# True parameter
true_p = 0.6
n_flips = 50

# Generate observation data
flips = np.random.binomial(1, true_p, n_flips)
n_heads = flips.sum()

print(f"Observation data: {n_flips} flips, {n_heads} heads")
print(f"MLE (Maximum Likelihood Estimate): p_hat = {n_heads/n_flips:.3f}")
print()

# Use Beta-Gamma relationship to sample from the posterior distribution
# Prior: Beta(2, 2)
# Posterior: Beta(2 + n_heads, 2 + n_flips - n_heads)

alpha_post = 2 + n_heads
beta_post = 2 + n_flips - n_heads

# Sample Beta distribution using Beta-Gamma relationship
def sample_beta(alpha, beta, n_samples=10000):
    """Sample from a Beta distribution using the Beta-Gamma relationship"""
    x = np.random.gamma(alpha, 1, n_samples)
    y = np.random.gamma(beta, 1, n_samples)
    return x / (x + y)

# Sample from the posterior
posterior_samples = sample_beta(alpha_post, beta_post, 10000)

# Posterior statistics
print("Posterior distribution statistics:")
print(f"  Posterior mean: {posterior_samples.mean():.4f}")
print(f"  Posterior standard deviation: {posterior_samples.std():.4f}")
print(f"  95% credible interval: [{np.percentile(posterior_samples, 2.5):.4f}, {np.percentile(posterior_samples, 97.5):.4f}]")

# Visualization
plt.figure(figsize=(10, 5))

plt.hist(posterior_samples, bins=50, density=True, alpha=0.7, 
         color='steelblue', edgecolor='black')

plt.axvline(true_p, color='green', linestyle='--', linewidth=2, label=f'True value p = {true_p}')
plt.axvline(n_heads/n_flips, color='red', linestyle=':', linewidth=2, label=f'MLE = {n_heads/n_flips:.3f}')
plt.axvline(posterior_samples.mean(), color='purple', linestyle='-.', linewidth=2, label=f'Posterior mean = {posterior_samples.mean():.3f}')

plt.xlabel('p')
plt.ylabel('Posterior Density')
plt.title(f'Bayesian Posterior Sampling (Beta({alpha_post}, {beta_post}))')
plt.legend()
plt.grid(alpha=0.3)

plt.tight_layout()
plt.show()
plt.close()

Click Run to execute code. Click the code area to edit.

Summary

This chapter provided hands-on programming practice for probability and statistical computation, starting from the most basic distribution sampling (uniform and normal distributions) and progressively delving into random seed control, random selection and shuffling operations, visualization of multiple probability distributions, descriptive statistics computation, covariance and correlation coefficient analysis, Monte Carlo simulation methods, and Bayesian posterior sampling techniques. These methods constitute the computational foundation of data analysis and machine learning. Uniform and normal distribution sampling are the starting points of random simulation; random seeds ensure experiment reproducibility; descriptive statistics provide a first-glance understanding of data characteristics; covariance and correlation coefficients quantify the strength of relationships between variables; Monte Carlo methods offer a general framework for approximate solutions to complex problems; and Bayesian posterior sampling achieves a complete description of uncertainty in parameter estimation.

Exercises

1. Use the Monte Carlo method to estimate $\int_0^1 x^2 dx$ and compare with the true value.

   Reference Answer

   import numpy as np
   
   n = 100000
   x = np.random.uniform(0, 1, n)
   estimate = np.mean(x**2)
   
   true_value = 1/3
   
   print(f"Monte Carlo estimate: {estimate:.6f}")
   print(f"True value: {true_value:.6f}")
   print(f"Error: {abs(estimate - true_value):.6f}")
   

   Click Run to execute code. Click the code area to edit.

2. Use the Monte Carlo method to verify the Central Limit Theorem (regardless of the original distribution, when the sample size is sufficiently large, the distribution of the sample mean approximately follows a normal distribution).

   Reference Answer

   import numpy as np
   import matplotlib.pyplot as plt
   
   # Sample from a uniform distribution
   n_samples = 10000
   sample_size = 30
   sample_means = []
   for _ in range(n_samples):
       sample = np.random.uniform(0, 1, sample_size)
       sample_means.append(np.mean(sample))
   sample_means = np.array(sample_means)
   
   # Theoretical values
   # X ~ U(0,1): E[X] = 0.5, Var[X] = 1/12
   # X̄: E[X̄] = 0.5, Var[X̄] = 1/(12*30)
   theoretical_mean = 0.5
   theoretical_std = np.sqrt(1 / (12 * sample_size))
   
   print(f"Mean of sample means: {sample_means.mean():.4f} (theoretical: {theoretical_mean})")
   print(f"Std of sample means: {sample_means.std():.4f} (theoretical: {theoretical_std:.4f})")
   
   # Visualization
   plt.hist(sample_means, bins=50, density=True, alpha=0.7)
   x = np.linspace(0.3, 0.7, 100)
   pdf = 1 / (theoretical_std * np.sqrt(2*np.pi)) * np.exp(-(x - theoretical_mean)**2 / (2 * theoretical_std**2))
   plt.plot(x, pdf, 'r-', linewidth=2, label='Theoretical Normal Distribution')
   plt.xlabel('Sample Mean')
   plt.ylabel('Density')
   plt.title('Central Limit Theorem Verification')
   plt.legend()
   plt.show()
   

   Click Run to execute code. Click the code area to edit.

Words: 3,762

Updated 2026-07-28