Chapter 3: Statistical Estimation and Hypothesis Testing

Introduction

In Chapter 1, we learned the fundamentals of descriptive statistics, which summarize the characteristics of the data at hand, and probability theory, which handles uncertainty mathematically. However, the data we can actually obtain is only a small part—a sample—of the entire group we are interested in (the population). We cannot inspect every product coming off a factory line, nor can we test a new drug on every patient in the world. This is where statistical inference becomes necessary: inferring the properties of a population from a limited sample.

This chapter covers the two pillars of statistical inference: estimation and hypothesis testing. In estimation, we learn how to infer population parameters (such as the mean and variance) from a sample. In hypothesis testing, we learn how to reach objective, data-driven conclusions to questions such as "does the new learning method really work?" or "is there a difference in quality between two production lines?" These are foundational techniques used throughout machine learning as well—in comparing model performance, running A/B tests, and evaluating the significance of features.

1. The Theory of Point Estimation and Interval Estimation

There are two main approaches to inferring a population parameter: estimating it with a single value, and estimating a range within which the value is likely to lie.

1.1 Point Estimation

Point estimation refers to inferring a single numerical value for a population parameter from sample data. The calculation method or statistic used for this inference is called an estimator, and the actual computed numerical value is called an estimate. For example, the sample mean $\bar{x}$ is an estimator of the population mean $\mu$, and an actually computed number such as "82.5 points" is the estimate.

Desirable Properties of Estimators

A good estimator is expected to have the following statistical properties.

• Unbiasedness: The property that the expected value of the estimator matches the true parameter. Mathematically, this is expressed as $E[\hat{\theta}] = \theta$. When an estimator is an unbiased estimator, repeating the same sampling procedure many times causes the average of the estimates to converge to the true value.
• Consistency: The property that as the sample size $n$ increases, the estimator converges in probability toward the true parameter.
• Efficiency: Among unbiased estimators, the property of having the smallest variance (the least spread in the estimates).

1.2 Interval Estimation

Point estimation tries to pin down a population parameter with a "single value," but since the estimate varies from sample to sample, it lacks information about how reliable that value is. Interval estimation addresses this by presenting a range in which the parameter is believed to lie, together with a certain level of confidence. A representative example of interval estimation is the confidence interval, which we cover in detail in Section 3. Interval estimation gives us a more informative inference than point estimation: "the true value is likely to be somewhere within this range."

1.3 Python Implementation: Confirming Point Estimation and Unbiasedness

Using adult male height data as an example, let's implement point estimation and use simulation to confirm that the unbiased variance converges toward the true population variance.

import numpy as np

np.random.seed(42)

# Population: adult male heights following a normal distribution with mean 170cm, std 8cm
population_mean = 170
population_std = 8

# Draw one sample and compute point estimates
sample = np.random.normal(population_mean, population_std, size=30)
point_estimate_mean = np.mean(sample)
point_estimate_var_biased = np.var(sample, ddof=0)      # Sample variance (biased)
point_estimate_var_unbiased = np.var(sample, ddof=1)     # Unbiased variance

print(f"Sample mean (point estimate of population mean): {point_estimate_mean:.3f} cm")
print(f"Sample variance (ddof=0, biased): {point_estimate_var_biased:.3f}")
print(f"Unbiased variance (ddof=1): {point_estimate_var_unbiased:.3f}")
print(f"True population mean: {population_mean} cm (normally unknown)")

print()
print("=== Confirming Unbiasedness of Estimators via Simulation ===")
n_simulations = 10000
sample_size = 30

biased_var_estimates = np.zeros(n_simulations)
unbiased_var_estimates = np.zeros(n_simulations)
mean_estimates = np.zeros(n_simulations)

for i in range(n_simulations):
    s = np.random.normal(population_mean, population_std, size=sample_size)
    mean_estimates[i] = np.mean(s)
    biased_var_estimates[i] = np.var(s, ddof=0)
    unbiased_var_estimates[i] = np.var(s, ddof=1)

true_variance = population_std ** 2
print(f"True population variance: {true_variance}")
print(f"Average of sample means ({n_simulations} trials): {np.mean(mean_estimates):.4f} (true value: {population_mean})")
print(f"Average of sample variance (ddof=0): {np.mean(biased_var_estimates):.4f} -> deviation from true value: {np.mean(biased_var_estimates) - true_variance:.4f}")
print(f"Average of unbiased variance (ddof=1): {np.mean(unbiased_var_estimates):.4f} -> deviation from true value: {np.mean(unbiased_var_estimates) - true_variance:.4f}")

Execution Result:

Sample mean (point estimate of population mean): 168.495 cm
Sample variance (ddof=0, biased): 50.113
Unbiased variance (ddof=1): 51.841
True population mean: 170 cm (normally unknown)

=== Confirming Unbiasedness of Estimators via Simulation ===
True population variance: 64
Average of sample means (10000 trials): 170.0013 (true value: 170)
Average of sample variance (ddof=0): 61.8140 -> deviation from true value: -2.1860
Average of unbiased variance (ddof=1): 63.9455 -> deviation from true value: -0.0545

2. The Principle and Implementation of Maximum Likelihood Estimation

Maximum Likelihood Estimation (MLE) is one of the most important methods for point estimation. It is based on the idea of "taking as the estimate of the true parameter the value that maximizes the probability (likelihood) of obtaining the observed data."

2.1 The Likelihood Function and Log-Likelihood

Suppose $n$ data points $x_1, x_2, \ldots, x_n$ are observed independently from a probability distribution with a parameter $\theta$. Viewing the probability (joint probability density) of obtaining this observed result as a function of $\theta$ gives us the likelihood function.

$$L(\theta) = \prod_{i=1}^{n} f(x_i; \theta)$$

Here, $f(x; \theta)$ is the probability density function with parameter $\theta$. Since working with the product directly is numerically inconvenient (it becomes an extremely small number and is prone to underflow), it is common to maximize the log-likelihood, obtained by taking the logarithm. Since the logarithm is a monotonically increasing function, the parameter that maximizes the log-likelihood is the same as the one that maximizes the likelihood function.

$$\ell(\theta) = \log L(\theta) = \sum_{i=1}^{n} \log f(x_i; \theta)$$

2.2 Finding the Maximum Likelihood Estimator

The Maximum Likelihood Estimator $\hat{\theta}_{MLE}$ is defined as the value of $\theta$ that maximizes the log-likelihood $\ell(\theta)$.

$$\hat{\theta}_{MLE} = \underset{\theta}{\arg\max}\ \ell(\theta)$$

For a normal distribution, differentiating the log-likelihood with respect to the mean $\mu$ and variance $\sigma^2$ and setting each to zero (solving the likelihood equations) yields the following closed-form solutions.

$$\hat{\mu}_{MLE} = \bar{x} = \frac{1}{n}\sum_{i=1}^{n}x_i \qquad \hat{\sigma}^2_{MLE} = \frac{1}{n}\sum_{i=1}^{n}(x_i - \bar{x})^2$$

2.3 Python Implementation: Maximum Likelihood Estimation via Numerical Optimization

Using the normal distribution and the exponential distribution as examples, let's numerically maximize the log-likelihood to estimate parameters, and confirm that the results match the closed-form solutions.

import numpy as np
from scipy import stats
from scipy.optimize import minimize

np.random.seed(42)

# --- Example 1: MLE for the mean and standard deviation of a normal distribution ---
true_mu, true_sigma = 170, 8
data = np.random.normal(true_mu, true_sigma, size=50)

def neg_log_likelihood_normal(params, data):
    mu, sigma = params
    if sigma <= 0:
        return np.inf
    return -np.sum(stats.norm.logpdf(data, loc=mu, scale=sigma))

initial_guess = [np.mean(data), np.std(data)]
result = minimize(neg_log_likelihood_normal, initial_guess, args=(data,), method='Nelder-Mead')
mle_mu, mle_sigma = result.x

# Closed-form solution (for a normal distribution, the MLE matches the sample mean and sample std (ddof=0))
closed_form_mu = np.mean(data)
closed_form_sigma = np.std(data, ddof=0)

# scipy.stats.norm.fit also performs the same maximum likelihood estimation internally
fit_mu, fit_sigma = stats.norm.fit(data)

print("=== Maximum Likelihood Estimation for Normal Distribution Parameters ===")
print(f"MLE via numerical optimization: mu={mle_mu:.4f}, sigma={mle_sigma:.4f}")
print(f"Closed-form solution:           mu={closed_form_mu:.4f}, sigma={closed_form_sigma:.4f}")
print(f"scipy.stats.norm.fit: mu={fit_mu:.4f}, sigma={fit_sigma:.4f}")

# --- Example 2: MLE for the rate parameter of an exponential distribution ---
print()
print("=== Maximum Likelihood Estimation for Exponential Distribution Parameter ===")
true_lambda = 0.5  # Mean waiting time = 1/lambda = 2 minutes
wait_times = np.random.exponential(scale=1/true_lambda, size=200)

def neg_log_likelihood_exp(params, data):
    lam = params[0]
    if lam <= 0:
        return np.inf
    return -np.sum(stats.expon.logpdf(data, scale=1/lam))

result_exp = minimize(neg_log_likelihood_exp, [1.0], args=(wait_times,), method='Nelder-Mead')
mle_lambda = result_exp.x[0]

# Closed-form solution for exponential distribution MLE: lambda_hat = 1 / sample mean
closed_form_lambda = 1 / np.mean(wait_times)

print(f"True parameter lambda: {true_lambda}")
print(f"MLE via numerical optimization: lambda={mle_lambda:.4f}")
print(f"Closed-form solution (1/sample mean): lambda={closed_form_lambda:.4f}")
print(f"Sample size: {len(wait_times)}, sample mean waiting time: {np.mean(wait_times):.4f} minutes")

Execution Result:

=== Maximum Likelihood Estimation for Normal Distribution Parameters ===
MLE via numerical optimization: mu=168.1962, sigma=7.3943
Closed-form solution:           mu=168.1962, sigma=7.3943
scipy.stats.norm.fit: mu=168.1962, sigma=7.3943

=== Maximum Likelihood Estimation for Exponential Distribution Parameter ===
True parameter lambda: 0.5
MLE via numerical optimization: lambda=0.4982
Closed-form solution (1/sample mean): lambda=0.4983
Sample size: 200, sample mean waiting time: 2.0069 minutes

All three calculation methods (numerical optimization, closed-form solution, and scipy's built-in function) give the same result. In the exponential distribution example as well, the MLE from numerical optimization closely matches the closed-form solution $\hat{\lambda} = 1/\bar{x}$, confirming that maximum likelihood estimation works exactly as the theory predicts.

3. Calculating and Interpreting Confidence Intervals

3.1 What Is a Confidence Interval?

The Confidence Interval (CI) is a representative technique for interval estimation, presenting a range in which the population parameter is expected to lie, together with a specified confidence level (often 95%). The confidence interval for the population mean $\mu$ is constructed around the sample mean $\bar{x}$, using the standard error $SE = s/\sqrt{n}$, which expresses the spread of the sample.

When the population standard deviation is unknown and the sample size is small (which is the case in most practical situations), the t-distribution is used instead of the standard normal distribution.

$$\bar{x} \pm t_{\alpha/2, n-1} \times \frac{s}{\sqrt{n}}$$

Here, $t_{\alpha/2, n-1}$ is the critical value from the t-distribution with $n-1$ degrees of freedom. As the sample size increases, the t-distribution approaches the standard normal distribution.

3.2 Python Implementation: Computing Confidence Intervals and Verifying Coverage

Using height data as an example, let's compute a 95% confidence interval, and then run 1,000 simulations to actually verify what the number "95%" means.

import numpy as np
from scipy import stats

np.random.seed(42)

# --- Confidence interval for the population mean (population std unknown, using the t-distribution) ---
population_mean = 170
population_std = 8
sample = np.random.normal(population_mean, population_std, size=30)

n = len(sample)
sample_mean = np.mean(sample)
sample_std = np.std(sample, ddof=1)  # Unbiased standard deviation
standard_error = sample_std / np.sqrt(n)

confidence_level = 0.95
degrees_of_freedom = n - 1

# Confidence interval using the t-distribution (scipy.stats.t.interval)
ci_lower, ci_upper = stats.t.interval(confidence_level, degrees_of_freedom,
                                       loc=sample_mean, scale=standard_error)

print("=== 95% Confidence Interval for the Population Mean ===")
print(f"Sample size: {n}")
print(f"Sample mean: {sample_mean:.3f} cm")
print(f"Sample standard deviation: {sample_std:.3f} cm")
print(f"Standard error: {standard_error:.3f} cm")
print(f"Degrees of freedom: {degrees_of_freedom}")
print(f"95% confidence interval: [{ci_lower:.3f}, {ci_upper:.3f}] cm")

# Verification by manual calculation
t_critical = stats.t.ppf(1 - (1 - confidence_level) / 2, degrees_of_freedom)
margin_of_error = t_critical * standard_error
print(f"\nVerification by manual calculation:")
print(f"t critical value (alpha=0.05, df={degrees_of_freedom}): {t_critical:.4f}")
print(f"Margin of error: {margin_of_error:.3f}")
print(f"Confidence interval: [{sample_mean - margin_of_error:.3f}, {sample_mean + margin_of_error:.3f}]")

# --- Verifying the coverage probability of confidence intervals via simulation ---
print()
print("=== Confirming the Meaning of the Confidence Interval's \"95%\" via Simulation ===")
n_simulations = 1000
contains_true_mean = 0

for _ in range(n_simulations):
    s = np.random.normal(population_mean, population_std, size=30)
    s_mean = np.mean(s)
    s_se = np.std(s, ddof=1) / np.sqrt(len(s))
    lo, hi = stats.t.interval(0.95, len(s) - 1, loc=s_mean, scale=s_se)
    if lo <= population_mean <= hi:
        contains_true_mean += 1

coverage_rate = contains_true_mean / n_simulations
print(f"Number of intervals (out of {n_simulations}) containing the true population mean ({population_mean}): {contains_true_mean}")
print(f"Coverage rate: {coverage_rate * 100:.2f}% (theoretical value: 95%)")

Execution Result:

=== 95% Confidence Interval for the Population Mean ===
Sample size: 30
Sample mean: 168.495 cm
Sample standard deviation: 7.200 cm
Standard error: 1.315 cm
Degrees of freedom: 29
95% confidence interval: [165.806, 171.183] cm

Verification by manual calculation:
t critical value (alpha=0.05, df=29): 2.0452
Margin of error: 2.689
Confidence interval: [165.806, 171.183]

=== Confirming the Meaning of the Confidence Interval's "95%" via Simulation ===
Number of intervals (out of 1000) containing the true population mean (170): 948
Coverage rate: 94.80% (theoretical value: 95%)

4. The Framework of Hypothesis Testing (Null Hypothesis, Alternative Hypothesis, p-value)

Hypothesis testing is a procedure for statistically judging the validity of a hypothesis about a population, based on sample data.

4.1 The Null Hypothesis and Alternative Hypothesis

• Null Hypothesis ($H_0$): A hypothesis such as "there is no difference" or "there is no effect"—the opposite of the claim we want to verify, and the hypothesis that is expected to be rejected.
• Alternative Hypothesis ($H_1$ or $H_a$): A hypothesis such as "there is a difference" or "there is an effect"—the claim we actually want to make.

In hypothesis testing, rather than directly proving that $H_1$ is true, we evaluate "how low is the probability of obtaining the observed data (or data even more extreme), assuming $H_0$ is true?" and reject $H_0$ if that probability is sufficiently low. This is a logical structure similar to proof by contradiction.

4.2 Significance Level and p-value

• Significance Level ($\alpha$): The threshold for the probability we are willing to accept of mistakenly rejecting a true null hypothesis. Conventionally, 0.05 (5%) or 0.01 (1%) is used.
• Test Statistic: A statistic computed from the sample data (such as a t-value or chi-squared value) used to judge the hypothesis.
• p-value: The probability, assuming the null hypothesis is true, of obtaining a test statistic as extreme as, or more extreme than, the one actually observed.

The decision rule is simple. If the p-value is below the significance level $\alpha$, we reject the null hypothesis and support the alternative hypothesis. If it is not below $\alpha$, we conclude that "there is insufficient evidence to reject the null hypothesis" (note that this does not mean "the null hypothesis has been proven true").

4.3 Type I and Type II Errors

A Type I Error is the mistake of judging "there is a difference" when there really is none (a false positive), and its probability is exactly the significance level $\alpha$. A Type II Error is the mistake of judging "there is no difference" when there really is one (a false negative). These two types of errors trade off against each other: tightening $\alpha$ tends to increase the Type II error rate.

4.4 Two-tailed Tests and One-tailed Tests

A Two-tailed Test verifies "there is a difference (regardless of direction)," with the alternative hypothesis taking the form $\mu \neq \mu_0$. A One-tailed Test verifies a specific direction of effect, with the alternative hypothesis taking the form $\mu > \mu_0$ or $\mu < \mu_0$. In most analyses, the two-tailed test is recommended as the default, to avoid arbitrary conclusions.

4.5 Python Implementation: Understanding the Hypothesis Testing Framework via Manual Calculation

Let's manually work through the procedure of a one-sample t-test to check whether the average score of students who received a new teaching method differs from the conventional average of 75 points.

import numpy as np
from scipy import stats

np.random.seed(42)

# Null hypothesis H0: the mean score of students using the new teaching method is the same as the conventional 75 points
# Alternative hypothesis H1: the mean score of students using the new teaching method differs from 75 points (two-tailed test)
population_mean_h0 = 75
sample_scores = np.array([78, 82, 75, 88, 79, 85, 91, 76, 83, 80,
                           77, 86, 89, 74, 81, 84, 90, 78, 82, 87])

n = len(sample_scores)
sample_mean = np.mean(sample_scores)
sample_std = np.std(sample_scores, ddof=1)
standard_error = sample_std / np.sqrt(n)

# Manually compute the t-statistic
t_statistic = (sample_mean - population_mean_h0) / standard_error
degrees_of_freedom = n - 1

# p-value for a two-tailed test
p_value = 2 * (1 - stats.t.cdf(abs(t_statistic), degrees_of_freedom))

alpha = 0.05

print("=== Hypothesis Testing Framework: One-sample t-test ===")
print(f"Null hypothesis H0: population mean = {population_mean_h0}")
print(f"Alternative hypothesis H1: population mean != {population_mean_h0} (two-tailed test)")
print(f"Significance level alpha: {alpha}")
print(f"Sample size: {n}, sample mean: {sample_mean:.3f}, sample std: {sample_std:.3f}")
print(f"Test statistic t: {t_statistic:.4f}")
print(f"Degrees of freedom: {degrees_of_freedom}")
print(f"p-value: {p_value:.2e}")

if p_value < alpha:
    print(f"Conclusion: since the p-value is below the significance level, we reject the null hypothesis")
else:
    print(f"Conclusion: since the p-value is not below the significance level, we do not reject the null hypothesis")

# Verification against scipy.stats.ttest_1samp
t_check, p_check = stats.ttest_1samp(sample_scores, population_mean_h0)
print(f"\nVerification with scipy.stats.ttest_1samp: t={t_check:.4f}, p={p_check:.2e}")

Execution Result:

=== Hypothesis Testing Framework: One-sample t-test ===
Null hypothesis H0: population mean = 75
Alternative hypothesis H1: population mean != 75 (two-tailed test)
Significance level alpha: 0.05
Sample size: 20, sample mean: 82.250, sample std: 5.149
Test statistic t: 6.2968
Degrees of freedom: 19
p-value: 4.81e-06
Conclusion: since the p-value is below the significance level, we reject the null hypothesis

Verification with scipy.stats.ttest_1samp: t=6.2968, p=4.81e-06

The manual calculation and scipy's built-in function ttest_1samp give exactly matching results. The p-value (about 0.0000048) is far below the significance level of 0.05, so we reject the null hypothesis that "the mean score with the new teaching method is unchanged from before," concluding that there is a significant difference.

5. Practical t-tests, Chi-squared Tests, and F-tests

Let's now organize and implement three testing methods that are frequently used in practice, grouped by purpose.

5.1 t-test

The t-test uses the t-distribution to compare means. You need to choose the appropriate variant depending on whether the data comes from two independent groups or is paired (the same subjects measured before and after).

import numpy as np
from scipy import stats

np.random.seed(42)

def judge(p, alpha=0.05):
    return "significant difference" if p < alpha else "no significant difference"

print("=" * 50)
print("1) Independent two-sample t-test (comparing means of two independent groups)")
print("=" * 50)
# Compare test scores between two teaching methods (Method A and Method B)
group_a = np.array([72, 75, 78, 80, 74, 77, 81, 79, 76, 73])
group_b = np.array([80, 85, 83, 88, 82, 86, 90, 84, 87, 81])

t_stat, p_value = stats.ttest_ind(group_a, group_b)
print(f"Group A mean: {np.mean(group_a):.2f}, Group B mean: {np.mean(group_b):.2f}")
print(f"t-statistic: {t_stat:.4f}")
print(f"p-value: {p_value:.3e}")
print(f"Judgment (alpha=0.05): {judge(p_value)}")

print()
print("=" * 50)
print("2) Paired two-sample t-test (paired t-test)")
print("=" * 50)
# Skill test scores before and after training (same subjects)
before = np.array([65, 70, 68, 72, 75, 69, 71, 74, 66, 73])
after = np.array([72, 78, 74, 80, 82, 75, 79, 81, 73, 80])

t_stat_paired, p_value_paired = stats.ttest_rel(before, after)
print(f"Mean before training: {np.mean(before):.2f}, mean after training: {np.mean(after):.2f}")
print(f"Mean difference: {np.mean(after - before):.2f}")
print(f"t-statistic: {t_stat_paired:.4f}")
print(f"p-value: {p_value_paired:.3e}")
print(f"Judgment (alpha=0.05): {judge(p_value_paired)}")

Execution Result:

==================================================
1) Independent two-sample t-test (comparing means of two independent groups)
==================================================
Group A mean: 76.50, Group B mean: 84.60
t-statistic: -5.8105
p-value: 1.665e-05
Judgment (alpha=0.05): significant difference

==================================================
2) Paired two-sample t-test (paired t-test)
==================================================
Mean before training: 70.30, mean after training: 77.40
Mean difference: 7.10
t-statistic: -30.4286
p-value: 2.189e-10
Judgment (alpha=0.05): significant difference

5.2 Chi-square Test

The chi-square test is used for categorical data (frequency data). It evaluates the discrepancy between observed and expected frequencies as a chi-squared statistic. Representative uses are the Goodness-of-fit Test (whether observed data follows a specific theoretical distribution) and the Test of Independence (whether two categorical variables are associated). The test of independence uses a contingency table, which arranges categories in rows and columns.

import numpy as np
from scipy import stats

def judge(p, alpha=0.05):
    return "significant difference" if p < alpha else "no significant difference"

print("=" * 50)
print("3) Chi-square Goodness-of-fit test")
print("=" * 50)
# Test whether a die is fair (equal probability for each face)
observed = np.array([18, 22, 16, 24, 20, 20])  # Observed frequency for each face (120 rolls total)
expected = np.array([20, 20, 20, 20, 20, 20])   # Expected frequency for a fair die

chi2_stat, p_value_chi2 = stats.chisquare(observed, expected)
print(f"Observed frequencies: {observed}")
print(f"Expected frequencies: {expected}")
print(f"Chi-squared statistic: {chi2_stat:.4f}")
print(f"p-value: {p_value_chi2:.4f}")
print(f"Judgment (alpha=0.05): {judge(p_value_chi2)}")

print()
print("=" * 50)
print("4) Chi-square test of independence")
print("=" * 50)
# Test whether there is an association between gender and product preference (A/B/C)
# Rows: gender (male, female), columns: preferred product (A, B, C)
contingency_table = np.array([
    [30, 20, 15],   # Male
    [15, 25, 25],   # Female
])
chi2_ind, p_ind, dof_ind, expected_ind = stats.chi2_contingency(contingency_table)
print(f"Contingency table:\n{contingency_table}")
print(f"Expected frequencies:\n{np.round(expected_ind, 2)}")
print(f"Chi-squared statistic: {chi2_ind:.4f}")
print(f"Degrees of freedom: {dof_ind}")
print(f"p-value: {p_ind:.4f}")
print(f"Judgment (alpha=0.05): {judge(p_ind)}")

Execution Result:

==================================================
3) Chi-square Goodness-of-fit test
==================================================
Observed frequencies: [18 22 16 24 20 20]
Expected frequencies: [20 20 20 20 20 20]
Chi-squared statistic: 2.0000
p-value: 0.8491
Judgment (alpha=0.05): no significant difference

==================================================
4) Chi-square test of independence
==================================================
Contingency table:
[[30 20 15]
 [15 25 25]]
Expected frequencies:
[[22.5 22.5 20. ]
 [22.5 22.5 20. ]]
Chi-squared statistic: 8.0556
Degrees of freedom: 2
p-value: 0.0178
Judgment (alpha=0.05): significant difference

In the goodness-of-fit test, the p-value of 0.8491 is far above 0.05, so there is no evidence to reject the null hypothesis that "this die is fair." In the test of independence, on the other hand, the p-value of 0.0178 is below 0.05, so we can conclude there is a statistically significant association between gender and product preference.

5.3 F-test

The F-test is a test that relies on the fact that the ratio of two variances (the F-statistic) follows the F-distribution. It is used both to directly compare the spread (variance) of two groups, and within Analysis of Variance (ANOVA), which simultaneously compares the means of three or more groups.

import numpy as np
from scipy import stats

def judge(p, alpha=0.05):
    return "significant difference" if p < alpha else "no significant difference"

print("=" * 50)
print("5) F-test (comparing the variance of two groups)")
print("=" * 50)
# Compare the variability (variance) of product weight between two production lines
line_a = np.array([100.2, 99.8, 100.5, 99.6, 100.1, 100.3, 99.9, 100.0, 99.7, 100.4])
line_b = np.array([100.1, 98.5, 101.3, 99.0, 100.8, 98.9, 101.5, 99.2, 100.6, 99.5])

var_a = np.var(line_a, ddof=1)
var_b = np.var(line_b, ddof=1)
f_stat = var_b / var_a  # Put the larger variance in the numerator
df1 = len(line_b) - 1
df2 = len(line_a) - 1
p_value_f = 2 * min(stats.f.cdf(f_stat, df1, df2), 1 - stats.f.cdf(f_stat, df1, df2))

print(f"Line A variance: {var_a:.5f}, Line B variance: {var_b:.5f}")
print(f"F-statistic (variance B / variance A): {f_stat:.4f}")
print(f"Degrees of freedom: ({df1}, {df2})")
print(f"p-value (two-tailed): {p_value_f:.4f}")
print(f"Judgment (alpha=0.05): {judge(p_value_f)}")

print()
print("=" * 50)
print("6) One-way ANOVA (application of the F-test: comparing means of three or more groups)")
print("=" * 50)
# Compare plant growth for three fertilizers (A, B, C)
fertilizer_a = np.array([12, 14, 13, 15, 12])
fertilizer_b = np.array([18, 17, 19, 16, 18])
fertilizer_c = np.array([22, 24, 21, 23, 25])

f_stat_anova, p_value_anova = stats.f_oneway(fertilizer_a, fertilizer_b, fertilizer_c)
print(f"Fertilizer A mean: {np.mean(fertilizer_a):.2f}")
print(f"Fertilizer B mean: {np.mean(fertilizer_b):.2f}")
print(f"Fertilizer C mean: {np.mean(fertilizer_c):.2f}")
print(f"F-statistic: {f_stat_anova:.4f}")
print(f"p-value: {p_value_anova:.3e}")
print(f"Judgment (alpha=0.05): {judge(p_value_anova)}")

Execution Result:

==================================================
5) F-test (comparing the variance of two groups)
==================================================
Line A variance: 0.09167, Line B variance: 1.14044
F-statistic (variance B / variance A): 12.4412
Degrees of freedom: (9, 9)
p-value (two-tailed): 0.0009
Judgment (alpha=0.05): significant difference

==================================================
6) One-way ANOVA (application of the F-test: comparing means of three or more groups)
==================================================
Fertilizer A mean: 13.20
Fertilizer B mean: 17.60
Fertilizer C mean: 23.00
F-statistic: 65.7091
p-value: 3.431e-07
Judgment (alpha=0.05): significant difference

6. The Multiple Comparisons Problem and Bonferroni Correction

6.1 What Is the Multiple Comparisons Problem?

When performing a single test at a significance level of $\alpha = 0.05$, the probability of mistakenly rejecting a true null hypothesis (a Type I error) is exactly 5%. However, when we repeat multiple tests on the same or similar data, the probability that at least one of them will "come out significant" by chance rises sharply. This is called the Multiple Comparisons Problem.

When $n$ independent tests are each performed at significance level $\alpha$, the probability that at least one comes out significant by chance (the Family-Wise Error Rate) is given by:

$$P(\text{at least one false positive}) = 1 - (1-\alpha)^n$$

For example, if you perform 20 tests at $\alpha=0.05$, this probability reaches about 64%. "If you run enough tests, something will always come out significant" is a statistical trap, and this frequently causes problems in feature selection for machine learning and multivariate comparisons in A/B testing as well.

6.2 Bonferroni Correction

The Bonferroni correction is the simplest and most conservative way to handle the multiple comparisons problem. When performing $n$ tests, it tightens the significance level used for each individual test from $\alpha$ to $\alpha/n$, keeping the family-wise error rate at or below the original $\alpha$.

$$\alpha_{\text{corrected}} = \frac{\alpha}{n}$$

An equivalent and commonly used alternative is to multiply each individual p-value by $n$ to get an adjusted p-value, and compare it against the original $\alpha$.

6.3 Python Implementation: Demonstrating the Multiple Comparisons Problem and Bonferroni Correction

Let's repeat 20 tests on data where the null hypothesis is true (there really is no difference), and confirm that false positives occur without correction, and that the Bonferroni correction correctly suppresses them.

import numpy as np
from scipy import stats

np.random.seed(3)

print("=== Demonstrating the Multiple Comparisons Problem: significant by chance even when H0 is true ===")
# Perform 20 independent tests (all on data where the null hypothesis H0: no mean difference is true)
n_tests = 20
alpha = 0.05
p_values = []

for i in range(n_tests):
    # Sample both groups from the same distribution (no mean difference)
    group_x = np.random.normal(50, 10, size=25)
    group_y = np.random.normal(50, 10, size=25)
    _, p = stats.ttest_ind(group_x, group_y)
    p_values.append(p)

p_values = np.array(p_values)
n_significant_uncorrected = np.sum(p_values < alpha)

print(f"Number of tests performed: {n_tests} (all with a true null hypothesis)")
print(f"Significance level alpha: {alpha}")
print(f"Number of tests significant (p<0.05) without correction: {n_significant_uncorrected}")
print(f"Theoretical probability of at least one false positive: 1 - (1-{alpha})^{n_tests} = {1 - (1-alpha)**n_tests:.4f}")

print()
print("=== Applying the Bonferroni Correction ===")
bonferroni_alpha = alpha / n_tests
n_significant_bonferroni = np.sum(p_values < bonferroni_alpha)

print(f"Corrected significance level: {alpha}/{n_tests} = {bonferroni_alpha:.5f}")
print(f"Number of tests significant after Bonferroni correction: {n_significant_bonferroni}")

print()
print("=== A Custom Function Equivalent to scipy/statsmodels' Bonferroni Correction ===")
def bonferroni_correction(p_values, alpha=0.05):
    """
    Apply the Bonferroni correction

    Parameters:
    -----------
    p_values : array-like
        List of p-values for each test
    alpha : float
        Significance level before correction

    Returns:
    --------
    corrected_alpha : float
        Significance level after correction
    rejected : ndarray of bool
        Whether each test is rejected (deemed significant)
    adjusted_p_values : ndarray
        Adjusted p-values (original p-value * number of tests, capped at 1.0)
    """
    p_values = np.asarray(p_values)
    n = len(p_values)
    corrected_alpha = alpha / n
    adjusted_p_values = np.minimum(p_values * n, 1.0)
    rejected = p_values < corrected_alpha
    return corrected_alpha, rejected, adjusted_p_values

corrected_alpha, rejected, adjusted_p = bonferroni_correction(p_values, alpha=0.05)
print(f"Corrected significance level: {corrected_alpha:.5f}")
print(f"Number of tests rejected: {np.sum(rejected)}")
print(f"Minimum adjusted p-value: {np.min(adjusted_p):.4f}")

Execution Result:

=== Demonstrating the Multiple Comparisons Problem: significant by chance even when H0 is true ===
Number of tests performed: 20 (all with a true null hypothesis)
Significance level alpha: 0.05
Number of tests significant (p<0.05) without correction: 2
Theoretical probability of at least one false positive: 1 - (1-0.05)^20 = 0.6415

=== Applying the Bonferroni Correction ===
Corrected significance level: 0.05/20 = 0.00250
Number of tests significant after Bonferroni correction: 0

=== A Custom Function Equivalent to scipy/statsmodels' Bonferroni Correction ===
Corrected significance level: 0.00250
Number of tests rejected: 0
Minimum adjusted p-value: 0.9014

7. Summary and Next Steps

In this chapter, we learned about estimation and hypothesis testing—the core of statistical inference for inferring the properties of a population from a limited sample.

Next Steps

Once you have mastered estimation and testing—the foundations of statistical inference—the next step is to move on to Bayesian statistics and the application of statistics to machine learning algorithms themselves. Continue with the following chapters in this series to put your statistical thinking to work in machine learning practice.

Practice Problems

import numpy as np
from scipy import stats

data = np.array([48.2, 51.5, 49.8, 50.3, 52.1, 47.9,
                  50.7, 49.5, 51.2, 48.8, 50.0, 51.8])

n = len(data)
mean = np.mean(data)
var = np.var(data, ddof=1)
std = np.std(data, ddof=1)
se = std / np.sqrt(n)

ci_lower, ci_upper = stats.t.interval(0.95, n - 1, loc=mean, scale=se)

print(f"Sample size: {n}")
print(f"Sample mean: {mean:.3f} mm")
print(f"Unbiased variance: {var:.3f}")
print(f"Unbiased standard deviation: {std:.3f} mm")
print(f"Standard error: {se:.3f} mm")
print(f"95% confidence interval: [{ci_lower:.3f}, {ci_upper:.3f}] mm")

import numpy as np
from scipy import stats

before_method = np.array([62, 65, 70, 68, 72, 66, 69, 71, 64, 67])
new_method = np.array([70, 75, 78, 74, 80, 76, 77, 79, 73, 75])

# Since we are comparing the means of two independent groups, use an independent two-sample t-test
t_stat, p_value = stats.ttest_ind(before_method, new_method)

print(f"Conventional method mean contracts: {np.mean(before_method):.2f}")
print(f"New method mean contracts: {np.mean(new_method):.2f}")
print(f"t-statistic: {t_stat:.4f}")
print(f"p-value: {p_value:.3e}")

alpha = 0.05
if p_value < alpha:
    print("Conclusion: there is a statistically significant difference")
else:
    print("Conclusion: there is no statistically significant difference")

import numpy as np

p_values = np.array([0.001, 0.012, 0.024, 0.038, 0.049])
alpha = 0.05
n_tests = len(p_values)

corrected_alpha = alpha / n_tests
adjusted_p_values = np.minimum(p_values * n_tests, 1.0)
rejected = p_values < corrected_alpha

print(f"Corrected significance level: {alpha}/{n_tests} = {corrected_alpha}")
for i, (p, adj, rej) in enumerate(zip(p_values, adjusted_p_values, rejected), start=1):
    print(f"Test {i}: p={p:.3f}, adjusted p-value={adj:.3f}, "
          f"judgment={'significant' if rej else 'not significant'}")