SciPy - Stats

SciPy - Statistics (scipy.stats)

The scipy.stats module is a comprehensive library for statistical functions and probability distributions. It contains a vast number of tools for statistical analysis, hypothesis testing, and working with random variables. You can get a detailed list of its capabilities using help(scipy.stats) in a Python console or by referring to the official SciPy documentation.

Probability distributions in scipy.stats are represented by objects. Each univariate distribution is an instance of a base class, categorized as continuous or discrete:

Class	Description
rv_continuous	A generic base class for continuous random variables, intended for subclassing to define specific distributions (e.g., normal, uniform, chi-squared).
rv_discrete	A generic base class for discrete random variables, intended for subclassing (e.g., binomial, poisson, geometric).
rv_histogram	Generates a distribution defined by a histogram of observed data.

Working with Continuous Random Variables: Normal Distribution

A continuous random variable can take any value within a given range. The normal (or Gaussian) distribution is one of the most common. In scipy.stats, it’s represented by norm. The loc keyword argument specifies the mean (μ), and scale specifies the standard deviation (σ).

The norm object (and other distribution objects) provides methods like:

• .cdf(x): Cumulative Distribution Function, $P(X \le x)$.
• .pdf(x): Probability Density Function (for continuous) or .pmf(x): Probability Mass Function (for discrete).
• .ppf(q): Percent Point Function (inverse of CDF), finds $x$ such that $P(X \le x) = q$.
• .rvs(size): Generates random variates (samples) from the distribution.
• .mean(), .std(), .var(): Calculate theoretical moments of the distribution.

Let’s compute the CDF for several points from a standard normal distribution (mean=0, std=1):

from scipy.stats import norm
import numpy as np

# Points at which to evaluate the CDF
points = np.array([-2., -1., 0., 1., 2., 3.])

# Calculate CDF values (default norm is standard normal: loc=0, scale=1)
cdf_values = norm.cdf(points)
print(f"CDF values for points {points}:\n{cdf_values}")

# Example for a normal distribution with mean=5, std=2
cdf_custom_norm = norm.cdf(points, loc=5, scale=2)
print(f"\nCDF for N(5,2) at points {points}:\n{cdf_custom_norm}")

The output will be:

CDF values for points [-2. -1.  0.  1.  2.  3.]:
[0.02275013 0.15865525 0.5        0.84134475 0.97724987 0.9986501 ]

CDF for N(5,2) at points [-2. -1.  0.  1.  2.  3.]:
[0.00023263 0.0013499  0.00620967 0.02275013 0.0668072  0.15865525]

To find the median (50th percentile) of a distribution, use the Percent Point Function (PPF):

from scipy.stats import norm

median_std_normal = norm.ppf(0.5) # For standard normal (loc=0, scale=1)
print(f"Median of standard normal distribution: {median_std_normal}")

median_custom_normal = norm.ppf(0.5, loc=10, scale=3)
print(f"Median of N(10,3) distribution: {median_custom_normal}")

Output:

Median of standard normal distribution: 0.0
Median of N(10,3) distribution: 10.0

To generate random samples from a distribution, use the rvs() method with the size argument. For reproducible results, set a seed using np.random.seed() or pass a random_state argument to rvs().

from scipy.stats import norm
import numpy as np

np.random.seed(42) # For reproducibility
random_samples = norm.rvs(loc=0, scale=1, size=5)
print(f"Random samples from N(0,1): {random_samples}")

Output (will be consistent due to the seed):

Random samples from N(0,1): [ 0.49671415 -0.1382643   0.64768854  1.52302986 -0.23415337]

Working with Discrete Random Variables: Binomial Distribution

The binomial distribution models the number of successes in a fixed number n of independent Bernoulli trials, each with a probability of success p. In scipy.stats, it’s binom.

from scipy.stats import binom
import numpy as np

# Parameters for binomial distribution
n_trials = 10  # Number of trials
p_success = 0.3 # Probability of success in each trial

# Probability Mass Function (PMF): P(X=k)
# Probability of exactly 5 successes in 10 trials
k_successes = 5
pmf_value = binom.pmf(k_successes, n_trials, p_success)
print(f"PMF P(X={k_successes} | n={n_trials}, p={p_success}): {pmf_value:.4f}")

# Cumulative Distribution Function (CDF): P(X <= k)
# Probability of 5 or fewer successes
cdf_value = binom.cdf(k_successes, n_trials, p_success)
print(f"CDF P(X<={k_successes} | n={n_trials}, p={p_success}): {cdf_value:.4f}")

# Generate random samples
np.random.seed(123)
random_binom_samples = binom.rvs(n_trials, p_success, size=8)
print(f"Random samples from Binom({n_trials}, {p_success}): {random_binom_samples}")

Output:

PMF P(X=5 | n=10, p=0.3): 0.1029
CDF P(X<=5 | n=10, p=0.3): 0.9527
Random samples from Binom(10, 0.3): [2 3 4 1 3 2 2 4]

Descriptive Statistics

scipy.stats also offers functions for descriptive statistics on datasets (typically NumPy arrays). Some key functions include:

Function	Description
describe(array, ...)	Computes several descriptive statistics of the input array (e.g., mean, variance, skewness, kurtosis).
gmean(array, ...)	Computes the geometric mean.
hmean(array, ...)	Computes the harmonic mean.
kurtosis(array, ...)	Computes the kurtosis (measure of tailedness). Fisher’s definition (normal ==> 0.0) is default.
mode(array, ...)	Returns the modal value (most frequent) and its count.
skew(array, ...)	Computes the skewness (measure of asymmetry).
iqr(array, ...)	Computes the Interquartile Range (IQR = 75th percentile - 25th percentile).
sem(array, ...)	Calculates the Standard Error of the Mean (SEM).
zscore(array, ...)	Calculates the z-score for each value relative to the sample mean and standard deviation.

Example using describe and other basic statistics:

from scipy import stats
import numpy as np

data = np.array([1, 2, 2, 3, 3, 3, 4, 4, 5, 10])

print(f"Data: {data}")
print(f"Mean: {np.mean(data):.2f}, Median: {np.median(data):.2f}")
print(f"Standard Deviation: {np.std(data):.2f}, Variance: {np.var(data):.2f}")

# Using scipy.stats.describe
descriptive_stats = stats.describe(data)
print("\nDescriptive statistics from stats.describe():")
print(f"  Number of observations: {descriptive_stats.nobs}")
print(f"  Min and max: {descriptive_stats.minmax}")
print(f"  Mean: {descriptive_stats.mean:.2f}")
print(f"  Variance: {descriptive_stats.variance:.2f}")
print(f"  Skewness: {descriptive_stats.skewness:.2f}")
print(f"  Kurtosis: {descriptive_stats.kurtosis:.2f}")

Output:

Data: [ 1  2  2  3  3  3  4  4  5 10]
Mean: 3.70, Median: 3.00
Standard Deviation: 2.26, Variance: 5.11

Descriptive statistics from stats.describe():
  Number of observations: 10
  Min and max: (1, 10)
  Mean: 3.70
  Variance: 5.68
  Skewness: 1.70
  Kurtosis: 2.64

Note: np.var calculates population variance by default, while stats.describe (and np.var(ddof=1)) calculates sample variance. This explains small differences in variance values if defaults are used.

Hypothesis Testing: T-tests

T-tests are used to determine if there is a significant difference between the means of two groups, or between a sample mean and a known population mean. It’s crucial to understand their assumptions (e.g., normality, independence of samples, equal variances for some tests).

ttest_1samp: One-Sample T-test

This test checks if the mean of a single sample is significantly different from a known or hypothesized population mean. It’s a two-sided test for the null hypothesis that the sample’s expected value (mean) equals the given population mean (popmean).

from scipy import stats
import numpy as np

# Generate some sample data (e.g., from a normal distribution)
np.random.seed(101)
sample_data = stats.norm.rvs(loc=5.5, scale=1.5, size=30)

# Test if the sample mean is significantly different from 5.0
population_mean_h0 = 5.0
t_statistic, p_value = stats.ttest_1samp(sample_data, population_mean_h0)

print(f"Sample mean: {np.mean(sample_data):.2f}")
print(f"T-statistic: {t_statistic:.3f}")
print(f"P-value: {p_value:.3f}")

alpha = 0.05 # Significance level
if p_value < alpha:
    print(f"Reject null hypothesis: Sample mean is significantly different from {population_mean_h0}.")
else:
    print(f"Fail to reject null hypothesis: No significant difference from {population_mean_h0}.")

Output (will vary slightly due to random sampling, but seed makes it reproducible):

Sample mean: 5.60
T-statistic: 2.166
P-value: 0.038
Reject null hypothesis: Sample mean is significantly different from 5.0.

ttest_ind: Independent Two-Sample T-test

This test compares the means of two independent samples. The null hypothesis is that the two independent samples have identical average (expected) values. By default, it assumes that the populations have identical variances (equal_var=True). If variances are unequal, set equal_var=False to perform Welch’s t-test.

from scipy import stats
import numpy as np

np.random.seed(42)
# Sample 1: e.g., scores from group A
sample1 = stats.norm.rvs(loc=5, scale=2, size=100)
# Sample 2: e.g., scores from group B (potentially different mean)
sample2 = stats.norm.rvs(loc=5.8, scale=2, size=100)

# Perform t-test assuming equal variances (default)
t_stat_eq_var, p_val_eq_var = stats.ttest_ind(sample1, sample2)
print(f"Independent T-test (equal variances assumed):")
print(f"  T-statistic: {t_stat_eq_var:.3f}, P-value: {p_val_eq_var:.3f}")

# Perform Welch's t-test (variances not assumed equal)
# Let's create another sample with different variance for illustration
sample3 = stats.norm.rvs(loc=5.8, scale=4, size=100) # Different scale
t_stat_welch, p_val_welch = stats.ttest_ind(sample1, sample3, equal_var=False)
print(f"\nWelch's T-test (unequal variances assumed for sample1 vs sample3):")
print(f"  T-statistic: {t_stat_welch:.3f}, P-value: {p_val_welch:.3f}")

Output (seed-dependent):

Independent T-test (equal variances assumed):
  T-statistic: -2.599, P-value: 0.010

Welch's T-test (unequal variances assumed for sample1 vs sample3):
  T-statistic: -1.872, P-value: 0.063

Interpreting P-values: A small p-value (typically < 0.05) suggests that the observed data is unlikely if the null hypothesis is true, leading to its rejection. A larger p-value means the data is consistent with the null hypothesis.

scipy.stats offers many more distributions (e.g., uniform, poisson, chi2, f) and statistical tests (e.g., ANOVA (f_oneway), Kolmogorov-Smirnov test (kstest), Shapiro-Wilk test for normality (shapiro)). Always consult the documentation to choose the appropriate test and understand its assumptions.