Chapter 33: Descriptive Statistics & Distributions

Measures of Center

The center of a distribution is the value around which the data tends to cluster. Three measures are in common use, each answering a slightly different question.

The arithmetic mean is the sum of all values divided by their count:

The mean minimizes the sum of squared deviations, making it the natural center for least-squares methods. Its weakness is sensitivity to outliers: a single extreme value can drag it far from where most data lives.

The median is the middle value when observations are sorted. For even $n$, it is the average of the two middle values. The median minimizes the sum of absolute deviations and is far more resistant to outliers.

The mode is the most frequently occurring value. For continuous data it is the peak of the density; for discrete data it is the most common category. Distributions can be unimodal, bimodal, or multimodal, each with different modeling implications.

Measures of Spread

The variance quantifies average squared deviation from the mean:

The $n-1$ denominator — Bessel's correction — makes $s^2$ an unbiased estimator of the population variance $\sigma^2$. Using $n$ instead underestimates because we are measuring deviations from the sample mean, which is itself estimated from the same data.

The standard deviation $s = \sqrt{s^2}$ restores the original units. The coefficient of variation $\text{CV} = s / \bar{x}$ (defined when $\bar{x} \neq 0$) is a dimensionless measure useful when comparing spread across variables with different scales.

The interquartile range IQR $= Q_3 - Q_1$ measures the spread of the middle 50% of data and is robust to outliers. The mean absolute deviation (MAD) and the related median absolute deviation $\text{MAD} = \text{median}(|x_i - \text{median}(x)|)$ are also robust alternatives.

import numpy as np
from scipy import stats

rng = np.random.default_rng(42)
x = rng.exponential(scale=2, size=500)   # right-skewed data

print(f"Mean:   {x.mean():.3f}")
print(f"Median: {np.median(x):.3f}")
print(f"Std:    {x.std(ddof=1):.3f}")
print(f"IQR:    {stats.iqr(x):.3f}")
print(f"MAD:    {stats.median_abs_deviation(x):.3f}")

For this right-skewed exponential sample you will see that the mean exceeds the median — a hallmark of positive skew.

Measures of Shape

Skewness measures asymmetry. The (Fisher) skewness of a distribution is:

where $\mu_3 = E[(X-\mu)^3]$ is the third central moment. Positive skew (long right tail) is common in income, wait times, and file sizes. Negative skew appears in test scores with a ceiling effect.

Kurtosis measures tail heaviness relative to a normal distribution. The excess kurtosis is $\gamma_2 = \mu_4/\sigma^4 - 3$. A positive excess kurtosis (leptokurtic) means heavier tails than normal — relevant for financial returns, where extreme events occur more often than a Gaussian model predicts. A negative excess kurtosis (platykurtic) means lighter tails.

print(f"Skewness:        {stats.skew(x):.3f}")
print(f"Excess Kurtosis: {stats.kurtosis(x):.3f}")  # Fisher definition, excess

Pitfalls

• Units matter for variance: if <!--MATHBLOCK15--> is in meters, <!--MATHBLOCK16--> is in meters-squared. Always prefer <!--MATHBLOCK17--> for communication.
• Mean vs. median storytelling: median household income is more informative than mean income because a billionaire shifts the mean dramatically.
• Sample size for shape statistics: skewness and kurtosis estimates are noisy for small <!--MATHBLOCK18--> — meaningful only for <!--MATHBLOCK19-->.

Quantiles and the Five-Number Summary

The $p$-th quantile (or $100p$-th percentile) $Q(p)$ is the value below which a fraction $p$ of the data falls. The quartiles $Q_1 = Q(0.25)$, $Q_2 = Q(0.5)$ (the median), and $Q_3 = Q(0.75)$ divide the dataset into four equal parts.

The five-number summary — minimum, $Q_1$, median, $Q_3$, maximum — gives a compact, robust snapshot of any distribution. It is the basis of the box plot (box-and-whisker plot), where the box spans the IQR, the whiskers extend to the farthest points within $1.5 \times \text{IQR}$ of the quartiles, and any points beyond the whiskers are plotted individually as suspected outliers.

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

rng = np.random.default_rng(7)
x = rng.lognormal(mean=1.5, sigma=0.8, size=300)

q1, med, q3 = np.percentile(x, [25, 50, 75])
iqr = q3 - q1
print(f"Min={x.min():.2f}  Q1={q1:.2f}  Med={med:.2f}  Q3={q3:.2f}  Max={x.max():.2f}")
print(f"IQR = {iqr:.2f}")

fig, ax = plt.subplots(figsize=(6, 3))
ax.boxplot(x, vert=False, patch_artist=True)
ax.set_xlabel("Value")
ax.set_title("Box plot of log-normal sample")
plt.tight_layout()
plt.savefig("boxplot_example.png", dpi=120)

Interpolation Methods for Quantiles

NumPy's np.percentile and R's quantile() both expose multiple interpolation methods (linear, nearest, midpoint, lower, upper, hazen, etc.). The default linear method is standard for most purposes. For small samples it can matter which you choose — be explicit.

Robust Measures of Center

The trimmed mean with trim fraction $\alpha$ discards the lowest and highest $\alpha$ fraction before averaging:

where $x_{(1)} \leq x_{(2)} \leq \cdots \leq x_{(n)}$ are the order statistics. Common choices are $\alpha = 0.05$ or $0.10$. The trimmed mean is more resistant to outliers than the mean while remaining more efficient than the median for near-normal data.

Winsorization replaces extreme values with the boundary quantile values rather than discarding them:

from scipy.stats import trim_mean, mstats

print(f"Mean:          {x.mean():.3f}")
print(f"Trimmed mean:  {trim_mean(x, 0.10):.3f}")

# Winsorize at 5th and 95th percentiles
xw = mstats.winsorize(x, limits=[0.05, 0.05])
print(f"Winsorized mean: {xw.mean():.3f}")

Outlier Detection

Outliers are observations that deviate markedly from the bulk of the data. Several detection rules exist:

• Tukey fences: outlier if <!--MATHBLOCK16--> or <!--MATHBLOCK17--> ("inner fences"). Extreme outliers use <!--MATHBLOCK18-->.
• Z-score rule: outlier if <!--MATHBLOCK19-->. Works best for near-normal data; unreliable for heavy-tailed distributions because outliers inflate <!--MATHBLOCK20--> and <!--MATHBLOCK21-->.
• Modified Z-score (Iglewicz-Hoaglin): uses median and MAD instead: <!--MATHBLOCK22-->; flag <!--MATHBLOCK23-->.

med = np.median(x)
mad = stats.median_abs_deviation(x)
mod_z = 0.6745 * (x - med) / mad
outliers = x[np.abs(mod_z) > 3.5]
print(f"{len(outliers)} outliers detected by modified Z-score")

Why Robustness Matters

In many real-world domains — sensor data with occasional spikes, financial returns with fat tails, genomic data with batch effects — the mean and variance are poor summaries. Robust statistics give you stable estimates that do not chase noise. The key insight: the median has a breakdown point of 50%, meaning up to half the sample can be replaced by arbitrary values without the median becoming unbounded. The mean has a breakdown point of 0%.

Pitfall: Removing Outliers Without Investigation

Outlier detection rules identify candidates, not errors. Before discarding any point, ask: Is it a data-entry error? A measurement artifact? Or a genuinely extreme but real observation? Removing real extreme values biases your analysis. Document every decision.

The Empirical Distribution Function

Given a sample $x_1, \ldots, x_n$, the empirical cumulative distribution function (ECDF) is:

It is a step function that jumps by $1/n$ at each observed value. The Glivenko-Cantelli theorem guarantees that $\hat{F}_n \to F$ uniformly as $n \to \infty$, making the ECDF a nonparametric, consistent estimator of the true CDF with no smoothing assumptions.

The ECDF is especially valuable because it is exact — every observation is represented without binning or smoothing. It is the natural tool for computing empirical quantiles and for comparing two samples visually.

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)
x = rng.normal(loc=5, scale=2, size=200)

xs = np.sort(x)
ys = np.arange(1, len(x) + 1) / len(x)

fig, ax = plt.subplots(figsize=(7, 4))
ax.step(xs, ys, where='post', color='steelblue', lw=2, label='ECDF')
ax.set_xlabel('Value')
ax.set_ylabel('$\\hat{F}_n(x)$')
ax.set_title('Empirical CDF of Normal(5, 2) sample')
ax.legend()
plt.tight_layout()

Histograms

A histogram bins the data into intervals and counts (or normalizes to density) how many observations fall in each bin. The bin width $h$ is the critical tuning parameter:

• Too narrow: high variance, noisy appearance.
• Too wide: high bias, important structure hidden.

Common automatic rules for selecting $h$:

• Sturges: <!--MATHBLOCK8--> bins. Too few bins for large <!--MATHBLOCK9-->, poor for skewed data.
• Scott: <!--MATHBLOCK10-->. Optimal for normal data in mean squared error sense.
• Freedman-Diaconis: <!--MATHBLOCK11-->. Robust version of Scott's rule, better for skewed data.

Matplotlib's default is 'auto' which picks the better of Sturges and Freedman-Diaconis. Always try multiple bin widths.

fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharey=False)
for ax, rule in zip(axes, ['auto', 5, 100]):
    ax.hist(x, bins=rule, density=True, color='steelblue', alpha=0.7, edgecolor='white')
    ax.set_title(f'bins={rule}')
    ax.set_xlabel('Value')
plt.tight_layout()

Kernel Density Estimation

KDE smooths the discrete mass at each observation using a kernel function $K$:

The most common choice is the Gaussian kernel $K(u) = \frac{1}{\sqrt{2\pi}}e^{-u^2/2}$. The bandwidth $h$ plays the same role as bin width: larger $h$ means more smoothing. Silverman's rule of thumb is $h = 0.9 \min(\hat{\sigma}, \text{IQR}/1.34) n^{-1/5}$, which is optimal for normal reference distributions.

from scipy.stats import gaussian_kde

kde = gaussian_kde(x, bw_method='silverman')
xgrid = np.linspace(x.min() - 1, x.max() + 1, 400)

fig, ax = plt.subplots(figsize=(7, 4))
ax.hist(x, bins='fd', density=True, color='steelblue', alpha=0.4, label='Histogram')
ax.plot(xgrid, kde(xgrid), 'k-', lw=2, label='KDE (Silverman)')
ax.set_xlabel('Value')
ax.set_ylabel('Density')
ax.legend()

KDE Pitfalls

• KDE leaks probability mass beyond the natural support of a variable (e.g., negative values for a positive-only variable). Use a log-transform or a boundary-corrected kernel for such data.
• For multimodal distributions, a global bandwidth may over-smooth one peak while under-smoothing another. Adaptive KDE methods address this.

ECDF vs. Histogram vs. KDE: When to Use Each

• Use the ECDF for precise quantile reading, two-sample comparison, and formal tests (Kolmogorov-Smirnov). It makes no smoothing assumptions.
• Use histograms for a quick look at shape, especially in presentations, when exact bin counts matter, or when the audience is non-technical.
• Use KDE for smooth density visualization, overlaying multiple groups, and publication-quality figures. Remember it requires a bandwidth choice.

For diagnostic purposes — checking normality, comparing samples — the ECDF is often the most honest tool.

The Bernoulli Distribution

The simplest discrete distribution models a single binary trial with success probability $p$:

Mean: $\mu = p$. Variance: $\sigma^2 = p(1-p)$. It is the building block for everything that follows.

The Binomial Distribution

If we run $n$ independent Bernoulli($p$) trials and count successes, we get $X \sim \text{Binomial}(n, p)$:

Mean: $np$. Variance: $np(1-p)$. The binomial is appropriate when:

• The number of trials <!--MATHBLOCK12--> is fixed in advance.
• Each trial is independent.
• The success probability <!--MATHBLOCK13--> is constant across trials.

Typical use cases: number of defective items in a batch of $n$; number of users who click a button given a sample of $n$ visitors; number of heads in $n$ coin flips.

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

n, p = 20, 0.3
k = np.arange(0, n+1)
pmf = binom.pmf(k, n, p)

fig, ax = plt.subplots(figsize=(7, 4))
ax.vlines(k, 0, pmf, colors='steelblue', lw=3)
ax.set_title(f'Binomial({n}, {p}) PMF')
ax.set_xlabel('k')
ax.set_ylabel('P(X = k)')

print(f"Mean:     {binom.mean(n, p):.2f}")
print(f"Variance: {binom.var(n, p):.2f}")
print(f"P(X >= 10): {1 - binom.cdf(9, n, p):.4f}")

Normal Approximation

When $np \geq 5$ and $n(1-p) \geq 5$, the binomial is well approximated by $\text{Normal}(np, np(1-p))$. Apply the continuity correction: $P(X \leq k) \approx \Phi\left(\frac{k + 0.5 - np}{\sqrt{np(1-p)}}\right)$.

The Poisson Distribution

The Poisson distribution models the number of events in a fixed interval of time or space when events occur at a constant average rate $\lambda$ and independently:

Mean: $\lambda$. Variance: $\lambda$. The equality of mean and variance is a defining characteristic — if you observe $\text{Var}(X) \gg E[X]$ (overdispersion), a negative binomial model is often more appropriate.

Typical use cases: call center arrivals per hour; mutations per genome per generation; photons hitting a detector per second; insurance claims per month.

from scipy.stats import poisson

lambda_ = 4
k = np.arange(0, 15)
pmf = poisson.pmf(k, lambda_)

fig, ax = plt.subplots(figsize=(7, 4))
ax.vlines(k, 0, pmf, colors='coral', lw=3)
ax.set_title(f'Poisson($\\lambda={lambda_}$) PMF')
ax.set_xlabel('k')
ax.set_ylabel('P(X = k)')

print(f"Mean=Var: {poisson.mean(lambda_):.2f}")
print(f"P(X = 0): {poisson.pmf(0, lambda_):.4f}")
print(f"P(X > 8): {1 - poisson.cdf(8, lambda_):.4f}")

Poisson as a Limit of the Binomial

When $n$ is large and $p$ is small with $\lambda = np$ fixed, $\text{Binomial}(n, p) \to \text{Poisson}(\lambda)$. This is why Poisson models rare events well: think of $n$ people in a city each independently having a rare car accident with probability $p$; the total count is approximately Poisson($np$).

The Negative Binomial

For overdispersed count data, the Negative Binomial $\text{NB}(r, p)$ models the number of failures before the $r$-th success. Parameterized by mean $\mu$ and dispersion $r$:

As $r \to \infty$ the variance approaches $\mu$ and we recover the Poisson. Negative binomial regression is the standard replacement for Poisson regression in the presence of overdispersion.

from scipy.stats import nbinom

# mean=4, dispersion r=2 -> Var = 4 + 16/2 = 12
r, mu = 2, 4
p_nb = r / (r + mu)
k = np.arange(0, 20)
plt.figure(figsize=(7,4))
plt.vlines(k, 0, nbinom.pmf(k, r, p_nb), colors='purple', lw=3, label='NegBin')
plt.vlines(k, 0, poisson.pmf(k, mu), colors='coral', lw=3, alpha=0.5, label='Poisson')
plt.legend(); plt.xlabel('k'); plt.title('Poisson vs Negative Binomial (same mean, larger variance)')

The Normal Distribution

The normal (Gaussian) distribution $X \sim \mathcal{N}(\mu, \sigma^2)$ has PDF:

The standard normal $Z = (X-\mu)/\sigma \sim \mathcal{N}(0, 1)$. The 68-95-99.7 rule: approximately 68%, 95%, and 99.7% of a normal distribution falls within 1, 2, and 3 standard deviations of the mean.

The normal distribution's dominance in statistics comes from the Central Limit Theorem: sums (and averages) of large numbers of independent, finite-variance random variables converge in distribution to normal, regardless of the original distribution. This is why measurement error, sample means, and many aggregate quantities are approximately normal.

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

x = np.linspace(-4, 4, 400)
fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# PDF
for mu, sig, label in [(0,1,'N(0,1)'), (0,2,'N(0,4)'), (2,1,'N(2,1)')]:
    axes[0].plot(x*sig + mu*0 + mu, norm.pdf(x*sig + mu*0 + mu, mu, sig), label=label)
axes[0].set_title('Normal PDFs'); axes[0].legend()

# Q-Q plot check
rng = np.random.default_rng(0)
data = rng.normal(3, 1.5, 300)
from scipy.stats import probplot
probplot(data, dist='norm', plot=axes[1])
axes[1].set_title('Normal Q-Q plot')
plt.tight_layout()

The Q-Q Plot

A quantile-quantile (Q-Q) plot compares the quantiles of your sample to the theoretical quantiles of a reference distribution. If the sample follows the reference distribution, the points fall on a 45-degree line. Systematic S-curves indicate heavier tails; a bow shape indicates skewness.

The Exponential Distribution

The exponential distribution $X \sim \text{Exp}(\lambda)$ models the waiting time between events in a Poisson process:

Mean: $1/\lambda$. Variance: $1/\lambda^2$. Its defining property is memorylessness: $P(X > s + t \mid X > s) = P(X > t)$. The time remaining until the next event does not depend on how long you have already waited — which is realistic for certain failure modes but unrealistic for aging systems.

Typical use cases: service times in queueing models; time to first failure of electronic components; inter-arrival times in a Poisson process.

The Gamma and Weibull Distributions

The Gamma distribution generalizes the exponential. $X \sim \text{Gamma}(\alpha, \beta)$ is the sum of $\alpha$ independent $\text{Exp}(1/\beta)$ variables:

Mean: $\alpha\beta$. Variance: $\alpha\beta^2$. When $\alpha = 1$ it reduces to exponential; when $\alpha = n/2, \beta = 2$ it gives the chi-squared distribution with $n$ degrees of freedom, central to statistical testing.

The Weibull distribution models failure times with non-constant hazard rates:

The shape parameter $k$ controls the hazard: $k < 1$ means decreasing hazard (infant mortality); $k = 1$ is exponential (constant hazard); $k > 1$ means increasing hazard (wear-out). It is the workhorse of reliability engineering.

Log-Normal and Heavy-Tailed Distributions

If $\ln X \sim \mathcal{N}(\mu, \sigma^2)$, then $X$ is log-normal. It arises naturally whenever a positive quantity is the product of many independent factors (by CLT on the log scale). Common examples: incomes, stock prices, particle sizes, file sizes, reaction times.

The Pareto and power-law distributions model truly heavy-tailed phenomena — where extreme values are far more common than a normal or exponential model predicts:

If $\alpha \leq 2$ the variance is infinite; if $\alpha \leq 1$ even the mean is infinite. Web page requests, earthquake magnitudes, and city sizes follow power laws.

from scipy.stats import lognorm, expon, gamma, weibull_min
import numpy as np
import matplotlib.pyplot as plt

xg = np.linspace(0.01, 10, 500)
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(xg, expon.pdf(xg, scale=2),       label='Exp(rate=0.5)')
ax.plot(xg, gamma.pdf(xg, a=2, scale=1),  label='Gamma(2,1)')
ax.plot(xg, weibull_min.pdf(xg, c=2, scale=2), label='Weibull(k=2)')
ax.plot(xg, lognorm.pdf(xg, s=0.8, scale=np.exp(1)), label='LogNormal')
ax.set_xlabel('x'); ax.set_ylabel('f(x)')
ax.set_title('Positive continuous distributions')
ax.legend(); plt.tight_layout()

Choosing a Distribution: A Practical Guide

• Bounded integers: Binomial, Hypergeometric
• Unbounded counts: Poisson, Negative Binomial
• Symmetric, approximately normal: Normal
• Positive, right-skewed, moderate: Gamma, Log-normal
• Failure times, reliability: Weibull, Exponential
• Extreme values / fat tails: Pareto, Generalized Extreme Value
• Proportions / probabilities: Beta

Always start with domain knowledge (what process generated the data?) before fitting. Then use Q-Q plots and goodness-of-fit tests to validate.

Maximum Likelihood Estimation

Given a parametric family $f(x; \theta)$ and data $x_1, \ldots, x_n$, maximum likelihood estimation (MLE) finds the parameter vector $\hat{\theta}$ that maximizes the likelihood:

Working in log-space converts the product to a sum and is numerically stable. SciPy implements MLE for all major distributions via .fit(data), which returns the shape parameters, location (loc), and scale.

from scipy.stats import gamma, norm, expon
import numpy as np

rng = np.random.default_rng(7)
data = rng.gamma(shape=2.5, scale=3.0, size=1000)

# MLE — fix loc=0 for a proper positive distribution
a, loc, scale = gamma.fit(data, floc=0)
print(f"Estimated shape={a:.3f} (true 2.5), scale={scale:.3f} (true 3.0)")

# Log-likelihood at MLE
log_lik = gamma.logpdf(data, a, loc, scale).sum()
print(f"Log-likelihood: {log_lik:.2f}")

Pitfall: loc Parameter

SciPy's dist.fit() estimates loc (shift) by default. For naturally positive data like durations, fix floc=0. Failing to do so can give nonsensical results — e.g., expon.fit(data) will set loc = min(data) rather than 0.

Method of Moments

An alternative to MLE: set theoretical moments equal to sample moments and solve for parameters. For the Gamma($\alpha$, $\beta$):

Method of moments is faster but generally less efficient (higher variance) than MLE. Use it as a starting point or when you need closed-form estimates.

Goodness-of-Fit Diagnostics

Q-Q Plots

The Q-Q plot is the primary visual diagnostic. Plot the empirical quantiles against the theoretical quantiles of the fitted distribution. Deviations indicate where the model fails:

• Heavy right tail: points curve upward above the line at the right end.
• Light tail: points curve downward below the line.
• Systematic S-curve: both tails depart — the distribution is heavier-tailed than the fit.
• Outliers: isolated extreme points far from the line.

import matplotlib.pyplot as plt
from scipy.stats import probplot

fig, axes = plt.subplots(1, 2, figsize=(11, 4))

# Good fit: Gamma vs Gamma-distributed data
probplot(data, dist=gamma(a, loc, scale), plot=axes[0])
axes[0].set_title('Q-Q: data vs fitted Gamma (good fit)')

# Bad fit: Normal vs Gamma-distributed data
probplot(data, dist='norm', plot=axes[1])
axes[1].set_title('Q-Q: data vs Normal (poor fit)')

plt.tight_layout()

Formal Tests

• Kolmogorov-Smirnov (KS) test: <!--MATHBLOCK9-->. Tests whether the sample came from a specific, fully specified <!--MATHBLOCK10-->. Critical limitation: when parameters are estimated from the same data, the standard KS critical values are anti-conservative (too liberal). Use Lilliefors correction for normality testing.
• Anderson-Darling (AD) test: weighted version of KS that gives more weight to the tails:

• Shapiro-Wilk test: specifically designed for normality. Highly recommended over KS/AD for testing normality; powerful even for small samples.

from scipy.stats import kstest, anderson, shapiro

# KS test against the fitted gamma
ks_stat, ks_p = kstest(data, lambda x: gamma.cdf(x, a, loc, scale))
print(f"KS test:  stat={ks_stat:.4f}, p={ks_p:.4f}")

# Anderson-Darling for normal
ad_result = anderson(data, dist='norm')
print(f"AD stat:  {ad_result.statistic:.4f}")
print(f"Critical values: {dict(zip(ad_result.significance_level, ad_result.critical_values))}")

# Shapiro-Wilk on a small normal sample
normal_sample = np.random.default_rng(0).normal(0,1,50)
stat, p_val = shapiro(normal_sample)
print(f"\nShapiro-Wilk on N(0,1) sample: stat={stat:.4f}, p={p_val:.4f}")

AIC and BIC for Model Selection

When comparing multiple candidate distributions, use information criteria:

where $k$ is the number of parameters and $\ell(\hat{\theta})$ is the maximized log-likelihood. Lower is better. BIC penalizes complexity more heavily for large $n$, favoring parsimonious models.

def aic(log_lik, k): return 2*k - 2*log_lik
def bic(log_lik, k, n): return k*np.log(n) - 2*log_lik

for dist_name, dist_obj, params, k in [
    ('Gamma',    gamma,  (a, loc, scale),                        2),
    ('Normal',   norm,   (data.mean(), data.std(ddof=1)),        2),
    ('Exponential', expon, (0, data.mean()),                     1),
]:
    ll = dist_obj.logpdf(data, *params).sum()
    print(f"{dist_name:12s}  AIC={aic(ll,k):.1f}  BIC={bic(ll,k,len(data)):.1f}")

A difference of more than 2 AIC units is generally considered meaningful; differences greater than 10 are strong evidence.