Probability distributions are mathematical functions that describe the likelihood of different possible outcomes in a random process. SciPy’s stats module provides useful tools for generating samples from these distributions and fitting distribution models to observed data. In this article, we will explore some of these techniques that SciPy offers. 

Getting Started with SciPy Stats

First, let’s import the necessary modules:

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

# Set random seed for reproducibility
np.random.seed(42)

Generating Samples from Different Distributions

SciPy makes it easy to generate random samples from a wide variety of distributions using the rvs (random variates) method:

# Generate 1000 samples from different distributions
normal_samples = stats.norm.rvs(loc=0, scale=1, size=1000)  # Normal distribution
uniform_samples = stats.uniform.rvs(loc=0, scale=10, size=1000)  # Uniform distribution
poisson_samples = stats.poisson.rvs(mu=5, size=1000)  # Poisson distribution

# Print summary statistics
print("Normal Distribution:")
print(f"Mean: {np.mean(normal_samples):.4f}")
print(f"Std Dev: {np.std(normal_samples):.4f}\n")

print("Uniform Distribution:")
print(f"Mean: {np.mean(uniform_samples):.4f}")
print(f"Range: [{np.min(uniform_samples):.4f}, {np.max(uniform_samples):.4f}]\n")

print("Poisson Distribution:")
print(f"Mean: {np.mean(poisson_samples):.4f}")
print(f"Variance: {np.var(poisson_samples):.4f}")

When we run this code, let’s examine what the output tells us about our generated samples:

Normal Distribution:
Mean: 0.0193
Std Dev: 0.9787

Uniform Distribution:
Mean: 5.0365
Range: [0.0322, 9.9941]

Poisson Distribution:
Mean: 4.9290
Variance: 4.3920

These results show how well SciPy generates samples that match the theoretical properties of each distribution:

• The normal distribution samples have a mean very close to 0 (0.0193) and standard deviation near 1 (0.9787), matching our specified parameters (μ=0, σ=1)
• The uniform distribution samples span almost exactly our desired range of 0 to 10, with a mean close to the theoretical 5
• The Poisson distribution samples have a mean of 4.929, very close to our specified λ=5, and its variance (4.392) is also close to λ as expected for a Poisson distribution

Let’s visualize these distributions to better understand their shapes:

# Create figure with three subplots
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 5))

# Plot Normal distribution
ax1.hist(normal_samples, bins=30, density=True, alpha=0.7, color='skyblue')
x = np.linspace(-4, 4, 100)
ax1.plot(x, stats.norm.pdf(x, loc=0, scale=1), 'r-', lw=2, label='PDF')
ax1.set_title('Normal Distribution\n(μ=0, σ=1)', fontsize=18)
ax1.legend(fontsize=15)

# Plot Uniform distribution
ax2.hist(uniform_samples, bins=30, density=True, alpha=0.7, color='lightgreen')
x = np.linspace(-1, 11, 100)
ax2.plot(x, stats.uniform.pdf(x, loc=0, scale=10), 'r-', lw=2, label='PDF')
ax2.set_title('Uniform Distribution\n(range: 0-10)', fontsize=18)
ax2.legend(fontsize=15)

# Plot Poisson distribution
ax3.hist(poisson_samples, bins=range(0, 15), density=True, alpha=0.7, color='salmon')
x = np.arange(0, 15)
ax3.plot(x, stats.poisson.pmf(x, mu=5), 'ro-', label='PMF')
ax3.set_title('Poisson Distribution\n(λ=5)', fontsize=18)
ax3.legend(fontsize=15)

plt.tight_layout()
plt.show()

The visualizations reveal the characteristic shapes of each distribution:

• The normal distribution shows the classic bell curve, symmetric around its mean
• The uniform distribution demonstrates equal probability across its range (the flat red line)
• The Poisson distribution, being discrete, shows the probability mass concentrated around its mean of 5

Fitting Distributions to Data

Now that we understand how these distributions look, let’s tackle a common real-world problem: determining which distribution best fits an unknown dataset. We will achieve this by creating a mystery data set (which actually comes from a gamma distribution). We will examine its properties and then later visualize to confirm our findings. 

# Generate mystery data (actually from a gamma distribution)
mystery_data = stats.gamma.rvs(a=2, loc=1.5, scale=2, size=1000)

# Try fitting different distributions
distributions = [
    ('Normal', stats.norm),
    ('Gamma', stats.gamma),
    ('Log-normal', stats.lognorm),
    ('Exponential', stats.expon)
]

# Fit each distribution and compute goodness of fit
results = []
for name, dist in distributions:
    # Fit distribution
    params = dist.fit(mystery_data)
    
    # Perform Kolmogorov-Smirnov test
    ks_stat, p_value = stats.kstest(mystery_data, dist.name, args=params)
    
    print(f"\n{name} distribution:")
    print(f"Parameters: {params}")
    print(f"KS statistic: {ks_stat:.4f}")
    print(f"P-value: {p_value:.4f}")

Let’s examine the fitting results:

Normal distribution:
Parameters: (5.51, 2.86)
KS statistic: 0.1089
P-value: 0.0000

Gamma distribution:
Parameters: (1.86, 1.54, 2.14)
KS statistic: 0.0201
P-value: 0.8041

Log-normal distribution:
Parameters: (0.60, 0.73, 4.03)
KS statistic: 0.0358
P-value: 0.1500

Exponential distribution:
Parameters: (1.57, 3.94)
KS statistic: 0.1246
P-value: 0.0000

Let’s interpret these results:

1. Goodness of Fit

The Kolmogorov-Smirnov (KS) test tells us how well each distribution fits our data:

• Lower KS statistics indicate better fits
• Higher p-values (> 0.05) suggest we cannot reject that the data follows that distribution
• The Gamma distribution has the best fit (lowest KS statistic: 0.0201, highest p-value: 0.8041)

2. Distribution Parameters

Each distribution is fitted with different parameters:

• Normal: mean=5.51, std=2.86
• Gamma: shape=1.86, location=1.54, scale=2.14
• Log-normal: s=0.60, location=0.73, scale=4.03
• Exponential: location=1.57, scale=3.94

Let’s go one step further to visualize how well each distribution fits our mystery data:

# Plot the mystery data and fitted distributions
plt.figure(figsize=(12, 6))

# Plot histogram of mystery data
plt.hist(mystery_data, bins=50, density=True, alpha=0.7, color='gray', label='Data')

# Plot fitted distributions
x = np.linspace(min(mystery_data), max(mystery_data), 100)
for name, dist in distributions:
    params = dist.fit(mystery_data)
    plt.plot(x, dist.pdf(x, *params), label=f'Fitted {name}', alpha=0.7)

plt.title('Mystery Data with Fitted Distributions', fontsize=20)
plt.xlabel('Value')
plt.ylabel('Density')
plt.legend(fontsize=16)
plt.grid(True, alpha=0.3)
plt.show()

The plot confirms our statistical findings:

• The Gamma distribution (orange line) follows our data’s shape most closely
• The Normal distribution (blue line) is too symmetric for our right-skewed data
• The Exponential distribution (red line) doesn’t capture the peak well
• The Log-normal distribution (green line) provides a reasonable fit but not as good as Gamma

Conclusion & Key Takeaways

Through this article, we’ve explored how SciPy provides some useful tools for working with probability distributions and fitting them to real-world data. From generating random samples to performing distribution fitting with statistical validation, we’ve seen how SciPy combines ease of use with robust statistical capabilities. Few key takeaways: 

1. Sample Generation: SciPy’s stats module provides efficient methods for generating samples that accurately match theoretical distributions, making it valuable for simulation and testing.
2. Complementary Analysis Tools: The combination of visual inspection and statistical tests offers a comprehensive approach to distribution analysis. While visualizations help understand the shape and characteristics of your data, statistical tests provide quantitative validation.
3. Distribution Fitting Best Practices:
   • Use p-values (> 0.05) and KS statistics to evaluate goodness of fit
   • Test multiple distribution types when the underlying distribution is unknown
   • Always validate statistical results with visual confirmation
4. Statistical Validation: The Kolmogorov-Smirnov test serves as a reliable metric for quantifying how well different distributions fit your data, helping make objective comparisons between potential fits.
5. Real-World Applications: Practical data often follows complex distributions like Gamma or Log-normal rather than simple Normal distributions, making SciPy’s diverse distribution options particularly valuable.

The combination of statistical testing and visual inspection, powered by SciPy’s comprehensive tools, provides data scientists with a robust framework for understanding and modeling probabilistic data distributions.