Hands-on Practical: Simulating and Plotting Distributions

Simulation and visualization offer a comprehensive understanding of probability distributions. While mathematical formulas provide the theoretical foundation, observing these distributions in action through simulation and visualization builds stronger intuition. This is especially relevant when evaluating if a particular distribution accurately fits observed data in a machine learning context.

We'll use Python, specifically the scipy.stats module, along with numpy for numerical operations and matplotlib or seaborn for plotting (often used with plotly for interactive web visualizations).

First, ensure you have the necessary libraries installed and import them:

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import seaborn as sns

# Optional: Set a style for seaborn plots
sns.set_style("whitegrid")

Simulating and Visualizing a Discrete Distribution: Binomial

Let's start with the Binomial distribution, which models the number of successes $k$ in a fixed number $n$ of independent Bernoulli trials, each with a success probability $p$.

Imagine we flip a biased coin ($p=0.6$) 20 times ($n=20$) and repeat this experiment many times (say, 1000). We can simulate the number of heads obtained in each experiment using scipy.stats.binom.rvs (random variates).

# Parameters for the Binomial distribution
n_binom = 20  # number of trials
p_binom = 0.6  # probability of success
size_binom = 1000 # number of simulations (experiments)

# Simulate random samples from the Binomial distribution
binomial_samples = stats.binom.rvs(n=n_binom, p=p_binom, size=size_binom)

# Print the first 10 simulated outcomes
print(f"First 10 simulated outcomes (number of successes out of {n_binom}):")
print(binomial_samples[:10])

Now, let's visualize the results of our simulation using a histogram and compare it to the theoretical probability mass function (PMF). The PMF gives the probability of obtaining exactly $k$ successes.

# Calculate the theoretical PMF
k_values_binom = np.arange(0, n_binom + 1)
pmf_binom = stats.binom.pmf(k=k_values_binom, n=n_binom, p=p_binom)

# Create the plot using Matplotlib (or adapt for Plotly)
plt.figure(figsize=(10, 6))

# Plot the histogram of simulated data
# Use density=True to normalize the histogram so it compares with the PMF
plt.hist(binomial_samples, bins=k_values_binom, density=True, alpha=0.6, color='#4dabf7', label='Simulated Data (Histogram)')

# Plot the theoretical PMF
# Use 'o-' to plot points connected by a line
plt.plot(k_values_binom, pmf_binom, 'o-', color='#f03e3e', label='Theoretical PMF')

plt.xlabel(f"Number of Successes (k) in {n_binom} Trials")
plt.ylabel("Probability / Density")
plt.title(f"Binomial Distribution (n={n_binom}, p={p_binom}) Simulation vs. Theory")
plt.legend()
plt.grid(True)
plt.show() 

Here is a representation of the plot for the web using Plotly:

Comparison of 1000 simulated Binomial experiments (n=20, p=0.6) with the theoretical probability mass function. The histogram closely follows the shape predicted by the PMF.

Notice how the histogram of our simulated results closely approximates the shape of the theoretical PMF. As the number of simulations (size_binom) increases, this approximation generally becomes better.

Simulating and Visualizing a Continuous Distribution: Normal

Next, let's work with the Normal (Gaussian) distribution, characterized by its mean $\mu$ and standard deviation $\sigma$ (or variance $\sigma^2$). It's arguably the most common distribution in statistics and machine learning, often used to model natural phenomena or errors.

Let's simulate data from a standard Normal distribution ($\mu=0, \sigma=1$) and compare the histogram of the samples to the theoretical probability density function (PDF).

# Parameters for the Normal distribution
mu_norm = 0     # mean
sigma_norm = 1  # standard deviation
size_norm = 1000 # number of samples

# Simulate random samples from the Normal distribution
normal_samples = stats.norm.rvs(loc=mu_norm, scale=sigma_norm, size=size_norm)

# Print the first 10 simulated values
print(f"First 10 simulated values from N({mu_norm}, {sigma_norm**2}):")
print(normal_samples[:10])

# Generate points for the theoretical PDF curve
x_values_norm = np.linspace(mu_norm - 4*sigma_norm, mu_norm + 4*sigma_norm, 200)
pdf_norm = stats.norm.pdf(x_values_norm, loc=mu_norm, scale=sigma_norm)

# Create the plot
plt.figure(figsize=(10, 6))

# Plot the histogram of simulated data
plt.hist(normal_samples, bins=30, density=True, alpha=0.6, color='#74c0fc', label='Simulated Data (Histogram)')

# Plot the theoretical PDF
plt.plot(x_values_norm, pdf_norm, color='#f76707', linewidth=2, label='Theoretical PDF')

plt.xlabel("Value")
plt.ylabel("Density")
plt.title(f"Normal Distribution (μ={mu_norm}, σ={sigma_norm}) Simulation vs. Theory")
plt.legend()
plt.grid(True)
plt.show() # Again, adapt for Plotly in web deployment

Here's the corresponding Plotly JSON representation:

Histogram of 1000 samples drawn from a standard Normal distribution (N(0, 1)), overlaid with the theoretical probability density function (PDF). The simulated data closely follows the characteristic bell shape.

Again, the histogram of the simulated data provides a good approximation of the bell-shaped curve defined by the theoretical PDF. This visualization helps confirm that our random samples behave as expected for a Normal distribution.

Simulating and Visualizing Another Distribution: Poisson

Let's do one more: the Poisson distribution. It models the probability of a given number of events occurring in a fixed interval of time or space, given the average rate of occurrence $\lambda$. Example: number of emails arriving per hour.

Assume emails arrive at an average rate of $\lambda=5$ per hour. We'll simulate the number of emails received in 1000 different hours.

# Parameters for the Poisson distribution
lambda_poisson = 5 # average rate (e.g., emails per hour)
size_poisson = 1000 # number of simulations (hours observed)

# Simulate random samples from the Poisson distribution
poisson_samples = stats.poisson.rvs(mu=lambda_poisson, size=size_poisson)

# Print the first 10 simulated counts
print(f"First 10 simulated counts (events per interval, λ={lambda_poisson}):")
print(poisson_samples[:10])

# Calculate the theoretical PMF
# Choose a reasonable upper limit for k, e.g., mean + few standard deviations
k_max_poisson = int(lambda_poisson + 4 * np.sqrt(lambda_poisson))
k_values_poisson = np.arange(0, k_max_poisson + 1)
pmf_poisson = stats.poisson.pmf(k=k_values_poisson, mu=lambda_poisson)

# Create the plot
plt.figure(figsize=(10, 6))

# Plot the histogram of simulated data
# Align bins correctly for discrete data
bins_poisson = np.arange(poisson_samples.min(), poisson_samples.max() + 2) - 0.5
plt.hist(poisson_samples, bins=bins_poisson, density=True, alpha=0.6, color='#69db7c', label='Simulated Data (Histogram)')

# Plot the theoretical PMF
plt.plot(k_values_poisson, pmf_poisson, 'o-', color='#be4bdb', label='Theoretical PMF')

plt.xlabel(f"Number of Events (k) per Interval")
plt.ylabel("Probability / Density")
plt.title(f"Poisson Distribution (λ={lambda_poisson}) Simulation vs. Theory")
plt.legend()
plt.grid(True)
plt.xticks(k_values_poisson) # Ensure ticks are at integer values
plt.xlim(bins_poisson.min(), bins_poisson.max())
plt.show()

Comparison of 1000 simulated Poisson experiments (lambda=5) with the theoretical probability mass function. The histogram mirrors the PMF, peaking around the average rate lambda.

This practical exercise demonstrates how simulation helps bridge the gap between abstract distribution formulas and concrete data. By generating samples and plotting them against theoretical functions using libraries like SciPy and Matplotlib/Plotly, you gain a better feel for the shape, spread, and characteristics of common probability distributions. This skill is valuable for exploratory data analysis, model selection, and understanding the probabilistic foundations of many machine learning algorithms.