Python SciPy Statistics: Distributions and Hypothesis Testing
python scipy statistics distributions and hypothesis testing: Practical guide to using scipy.stats for probability distributions, distribution fitting, and hypothesis...
The Statistical Toolkit in scipy.stats
When a Python project needs to model a probability distribution or validate a statistical claim, the scipy.stats module is usually the first place to turn. It combines a large set of continuous and discrete distributions with a collection of hypothesis tests that cover the most common data science and engineering workflows. This article walks through python scipy statistics distributions and hypothesis testing in the order a working developer would actually use them: from generating and inspecting distributions, to fitting them on real data, then running the right test for a given comparison.
Every distribution in scipy.stats exposes a consistent set of methods: pdf for the density, cdf for the cumulative probability, ppf for the inverse CDF, and rvs for random sampling. The parameters vary by family, but loc and scale are shared by all continuous distributions. loc shifts the distribution along the real line and scale stretches it, which means a standard normal is stats.norm(loc=0, scale=1) and a normal with mean 5 and standard deviation 2 is stats.norm(loc=5, scale=2).
from scipy import stats normal = stats.norm(loc=0.0, scale=1.0) print(normal.pdf(0.0)) # density at 0 print(normal.cdf(1.96)) # P(X <= 1.96) print(normal.ppf(0.975)) # quantile for 0.975
The ppf method is the inverse of cdf, so it answers the question "which value leaves a given probability mass to the left?" That is the operation behind confidence intervals and critical values.
Generating Random Samples from Distributions
The rvs method produces random draws from a distribution. It accepts size for the number of samples and random_state to make the draw reproducible. Passing an integer to random_state seeds the underlying generator, which matters when a test or a simulation must be rerun with identical data.
samples = stats.norm.rvs( loc=10.0, scale=2.0, size=1000, random_state=42, )
For discrete families, the same pattern applies. A Poisson draw uses stats.poisson.rvs(mu=3.0, size=500, random_state=0), and a binomial draw uses stats.binom.rvs(n=10, p=0.4, size=500, random_state=0). The returned array is a plain NumPy array, so it can be passed directly to plotting, aggregation, or further statistical routines.
One common mistake is to reuse a single random_state value across several independent draws that are meant to be independent. Reusing the same seed produces identical sequences, which silently correlates the samples. Use different seeds, or pass None and let the generator choose a fresh state.
Fitting a Distribution to Observed Data
When the underlying distribution is unknown, the fit method estimates its parameters by maximum likelihood. For a normal family, stats.norm.fit(data) returns a tuple (loc, scale). For a gamma family it returns (shape, loc, scale), and the order of the returned parameters follows the order of the distribution's own parameter list.
observed = stats.norm.rvs(loc=5.0, scale=1.5, size=500, random_state=7) loc_hat, scale_hat = stats.norm.fit(observed) print(loc_hat, scale_hat)
The quality of the fit depends on the sample size and on whether the chosen family actually matches the data. fit will always return parameters, even when the family is a poor match, so it is worth checking the result against a goodness-of-fit test rather than trusting the fitted curve visually. For heavy-tailed data, a normal fit will systematically underestimate the tails; a Student's t or a generalized extreme value family may be more appropriate.
One-Sample Hypothesis Tests
The most common one-sample question is whether the mean of a sample differs from a known reference value. ttest_1samp answers that question under the assumption that the data are approximately normal.
measurements = [12.1, 11.9, 12.4, 12.0, 12.2, 11.8] result = stats.ttest_1samp(measurements, popmean=12.0) print(result.statistic) print(result.pvalue)
The result object carries statistic and pvalue. A small p-value indicates that the observed mean is unlikely under the null hypothesis that the population mean equals popmean. The test does not tell you how large the difference is; report the sample mean and a confidence interval separately.
When the question is about the shape of the distribution rather than its mean, kstest compares the empirical CDF of the sample against a reference CDF. For example, stats.kstest(data, "norm", args=(loc, scale)) tests whether the data follow a normal distribution with the given parameters. The same function works with any distribution name that scipy.stats knows, and it is the standard tool for a goodness-of-fit check after fit.
For a quick normality check, shapiro tests the null hypothesis that the data came from a normal distribution. It is more sensitive than kstest for small samples, but it only applies to normality, not to other families.
Comparing Two Samples
Two-sample comparisons appear constantly in engineering work: a control group versus a treatment group, measurements before and after a change, or two implementations producing latency samples. The choice of test depends on whether the samples are paired and whether the normality assumption is acceptable.
For independent samples with approximately normal distributions, ttest_ind compares the means. The equal_var parameter controls whether the test assumes equal variances; when in doubt, set it to False to use Welch's variant, which does not require equal variances.
group_a = stats.norm.rvs(loc=10.0, scale=2.0, size=50, random_state=1) group_b = stats.norm.rvs(loc=11.0, scale=2.0, size=50, random_state=2) t_result = stats.ttest_ind(group_a, group_b, equal_var=False) print(t_result.pvalue)
For paired data, such as the same machine measured before and after a configuration change, ttest_rel is the correct choice because it accounts for the correlation between the paired observations.
When the normality assumption is not safe, the Mann-Whitney U test (mannwhitneyu) compares the distributions of two independent samples without assuming normality. It tests whether one sample tends to produce larger values than the other. The alternative parameter accepts "two-sided", "less", or "greater", matching the direction of the hypothesis. The Wilcoxon signed-rank test (wilcoxon) is the paired counterpart.
mw_result = stats.mannwhitneyu( group_a, group_b, alternative="two-sided", ) print(mw_result.pvalue)
ANOVA and Categorical Tests
When more than two groups must be compared, running pairwise t-tests inflates the chance of a false positive. A one-way ANOVA, implemented as f_oneway, tests whether the group means are all equal in a single procedure.
group_c = stats.norm.rvs(loc=9.0, scale=2.0, size=50, random_state=3) anova = stats.f_oneway(group_a, group_b, group_c) print(anova.statistic) print(anova.pvalue)
ANOVA assumes independence, approximate normality within each group, and roughly equal variances across groups. If those assumptions are violated, a Kruskal-Wallis test (kruskal) provides a non-parametric alternative that compares the distributions of the groups.
For categorical data, chi2_contingency tests whether two categorical variables are independent. It takes a contingency table as a 2D array and returns the chi-square statistic, the p-value, the degrees of freedom, and the expected counts under independence.
observed = [[120, 80], [90, 110]] chi2, p, dof, expected = stats.chi2_contingency(observed) print(chi2, p, dof) print(expected)
The expected counts are useful for diagnosing which cells deviate most from independence. A low expected count in any cell makes the chi-square approximation unreliable; a common rule of thumb is that all expected counts should be at least 5, though the test still runs and returns a result regardless.
Assumptions, Multiple Comparisons, and Runtime Cost
Every test in scipy.stats makes assumptions, and the p-value is only meaningful when those assumptions hold. The t-tests assume normality, though they are reasonably robust to moderate departures when the sample size is large. The Mann-Whitney and Wilcoxon tests assume that the two distributions have the same shape under the null hypothesis; they are not a test of medians when the shapes differ. The chi-square test assumes independent observations and sufficient expected counts in each cell.
When several tests are run on the same data, the probability of at least one false positive grows with the number of tests. A Bonferroni correction, which multiplies each p-value by the number of comparisons, is the simplest safeguard, but it is conservative. For exploratory work, treat small p-values from many tests as candidates for confirmation on a fresh dataset rather than as proof.
The runtime cost of these routines is generally small for typical sample sizes. Distribution fitting is the most expensive operation because it solves a numerical optimization problem; for very large arrays, fitting a gamma or a beta distribution takes noticeably longer than fitting a normal. Hypothesis tests are dominated by the cost of sorting or summing the data, so they scale roughly linearly with sample size. When the same distribution object is used repeatedly, construct it once and reuse it instead of recreating it inside a loop, since the object caches internal parameterization work.
For reproducible results, always pass random_state to rvs and record the seed in the same place as the analysis code. That makes a simulation rerunnable and lets a reviewer regenerate the exact dataset that produced a reported p-value.