F Distribution in Python: Complete Guide

Key Insights

The F distribution is fundamental to variance comparison and ANOVA, arising naturally when you divide two chi-squared distributions by their respective degrees of freedom—master it and you unlock powerful hypothesis testing capabilities.
SciPy’s scipy.stats.f module provides everything you need: PDF, CDF, inverse functions, random sampling, and moment calculations, all accessible through a consistent, intuitive API.
Understanding the relationship between degrees of freedom and distribution shape is critical—higher degrees of freedom produce distributions that are more symmetric and concentrated, directly impacting your statistical inference.

Introduction to the F Distribution

The F distribution, named after Ronald Fisher, is a continuous probability distribution that emerges when you take the ratio of two independent chi-squared random variables, each divided by their degrees of freedom. Mathematically, if X₁ ~ χ²(d₁) and X₂ ~ χ²(d₂) are independent, then:

F = (X₁/d₁) / (X₂/d₂)

follows an F distribution with d₁ (numerator) and d₂ (denominator) degrees of freedom.

This distribution appears constantly in statistical practice. ANOVA uses it to test whether group means differ significantly. Regression analysis relies on it to assess overall model significance. When you need to compare variances between two populations, the F-test is your tool. Understanding how to work with this distribution in Python isn’t optional—it’s essential for anyone doing serious statistical work.

The two degrees of freedom parameters (often called dfn and dfd in Python) control the distribution’s shape. Low degrees of freedom produce heavily right-skewed distributions. As both parameters increase, the distribution becomes more symmetric and approaches normality.

Working with F Distribution in SciPy

SciPy’s scipy.stats.f module is your primary interface for F distribution calculations. It follows the same pattern as other continuous distributions in SciPy, making it easy to learn if you’re familiar with the library.

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

# Create an F distribution with dfn=5, dfd=20
dfn, dfd = 5, 20
f_dist = stats.f(dfn, dfd)

# Explore basic properties
print(f"Mean: {f_dist.mean():.4f}")
print(f"Variance: {f_dist.var():.4f}")
print(f"Standard Deviation: {f_dist.std():.4f}")

# Theoretical mean of F distribution: dfd / (dfd - 2) for dfd > 2
theoretical_mean = dfd / (dfd - 2)
print(f"Theoretical mean: {theoretical_mean:.4f}")

# Get support (range where PDF > 0)
print(f"Support: [{f_dist.support()[0]:.4f}, {f_dist.support()[1]}]")

The F distribution is always positive (support starts at 0) and extends to infinity. The mean exists only when dfd > 2, and the variance exists only when dfd > 4. Keep these constraints in mind when working with small sample sizes.

Probability Density and Cumulative Distribution Functions

The PDF tells you the relative likelihood of observing a particular F value. The CDF gives you the probability of observing a value less than or equal to a given point—essential for calculating p-values.

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

# Define different degree of freedom combinations
df_combinations = [(2, 10), (5, 20), (10, 50), (50, 100)]
x = np.linspace(0.01, 5, 500)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Plot PDFs
for dfn, dfd in df_combinations:
    pdf_values = stats.f.pdf(x, dfn, dfd)
    axes[0].plot(x, pdf_values, label=f'dfn={dfn}, dfd={dfd}', linewidth=2)

axes[0].set_xlabel('x')
axes[0].set_ylabel('Probability Density')
axes[0].set_title('F Distribution PDF')
axes[0].legend()
axes[0].set_ylim(0, 1)
axes[0].grid(True, alpha=0.3)

# Plot CDFs
for dfn, dfd in df_combinations:
    cdf_values = stats.f.cdf(x, dfn, dfd)
    axes[1].plot(x, cdf_values, label=f'dfn={dfn}, dfd={dfd}', linewidth=2)

axes[1].set_xlabel('x')
axes[1].set_ylabel('Cumulative Probability')
axes[1].set_title('F Distribution CDF')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Practical calculations
dfn, dfd = 5, 20
f_value = 2.5

# Probability of F <= 2.5
prob_less = stats.f.cdf(f_value, dfn, dfd)
print(f"P(F <= {f_value}): {prob_less:.4f}")

# Survival function: P(F > 2.5) - useful for one-tailed tests
prob_greater = stats.f.sf(f_value, dfn, dfd)
print(f"P(F > {f_value}): {prob_greater:.4f}")

# Percent point function (inverse CDF): find F value for given probability
critical_95 = stats.f.ppf(0.95, dfn, dfd)
print(f"95th percentile (critical value): {critical_95:.4f}")

The survival function (sf) is particularly useful because F-tests are typically one-tailed—you’re looking at the right tail to determine if your F statistic is unusually large.

Generating Random Samples

When running simulations or bootstrapping, you’ll need to generate random samples from F distributions. Both SciPy and NumPy provide this capability.

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

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

dfn, dfd = 5, 20
n_samples = 10000

# Method 1: SciPy
samples_scipy = stats.f.rvs(dfn, dfd, size=n_samples)

# Method 2: NumPy (often faster for large samples)
rng = np.random.default_rng(42)
samples_numpy = rng.f(dfn, dfd, size=n_samples)

# Compare sample statistics to theoretical values
print("SciPy samples:")
print(f"  Sample mean: {samples_scipy.mean():.4f}, Theoretical: {dfd/(dfd-2):.4f}")
print(f"  Sample variance: {samples_scipy.var():.4f}")

# Visualize: histogram vs theoretical PDF
fig, ax = plt.subplots(figsize=(10, 6))

# Histogram of samples
ax.hist(samples_scipy, bins=100, density=True, alpha=0.7, 
        label='Sampled data', color='steelblue', edgecolor='white')

# Theoretical PDF
x = np.linspace(0, 8, 500)
pdf_theoretical = stats.f.pdf(x, dfn, dfd)
ax.plot(x, pdf_theoretical, 'r-', linewidth=2.5, label='Theoretical PDF')

ax.set_xlabel('F value')
ax.set_ylabel('Density')
ax.set_title(f'F Distribution Samples vs Theoretical (dfn={dfn}, dfd={dfd})')
ax.legend()
ax.set_xlim(0, 8)
ax.grid(True, alpha=0.3)
plt.show()

For most applications, NumPy’s random generator is faster, especially for large sample sizes. Use SciPy when you need access to other distribution methods in the same workflow.

Statistical Inference with F Distribution

The F-test compares variances from two populations. Under the null hypothesis that both populations have equal variance, the ratio of sample variances follows an F distribution.

import numpy as np
from scipy import stats

np.random.seed(123)

# Generate two samples with different variances
sample1 = np.random.normal(loc=50, scale=10, size=30)  # variance = 100
sample2 = np.random.normal(loc=50, scale=15, size=25)  # variance = 225

# Calculate sample variances
var1 = np.var(sample1, ddof=1)
var2 = np.var(sample2, ddof=1)

print(f"Sample 1 variance: {var1:.4f}")
print(f"Sample 2 variance: {var2:.4f}")

# F-test: always put larger variance in numerator for two-tailed test
if var1 >= var2:
    f_statistic = var1 / var2
    dfn = len(sample1) - 1
    dfd = len(sample2) - 1
else:
    f_statistic = var2 / var1
    dfn = len(sample2) - 1
    dfd = len(sample1) - 1

print(f"\nF-statistic: {f_statistic:.4f}")
print(f"Degrees of freedom: dfn={dfn}, dfd={dfd}")

# Calculate p-value (two-tailed)
p_value = 2 * stats.f.sf(f_statistic, dfn, dfd)
print(f"Two-tailed p-value: {p_value:.4f}")

# Critical values for alpha = 0.05
alpha = 0.05
critical_lower = stats.f.ppf(alpha/2, dfn, dfd)
critical_upper = stats.f.ppf(1 - alpha/2, dfn, dfd)
print(f"\nCritical values (α=0.05): [{critical_lower:.4f}, {critical_upper:.4f}]")

# Decision
if p_value < alpha:
    print("Reject null hypothesis: variances are significantly different")
else:
    print("Fail to reject null hypothesis: no significant difference in variances")

# Confidence interval for variance ratio
ci_lower = f_statistic / stats.f.ppf(1 - alpha/2, dfn, dfd)
ci_upper = f_statistic / stats.f.ppf(alpha/2, dfn, dfd)
print(f"\n95% CI for variance ratio: [{ci_lower:.4f}, {ci_upper:.4f}]")

This manual implementation gives you full control. In practice, you might use scipy.stats.levene or scipy.stats.bartlett for variance homogeneity tests, but understanding the underlying F-test mechanics is crucial.

Practical Application: ANOVA F-Statistic

One-way ANOVA tests whether the means of three or more groups differ significantly. The F-statistic measures the ratio of between-group variance to within-group variance.

import numpy as np
from scipy import stats

np.random.seed(456)

# Three treatment groups with different means
group_a = np.random.normal(loc=20, scale=5, size=25)
group_b = np.random.normal(loc=23, scale=5, size=30)
group_c = np.random.normal(loc=26, scale=5, size=28)

# Manual ANOVA calculation
all_data = np.concatenate([group_a, group_b, group_c])
grand_mean = all_data.mean()
n_total = len(all_data)
k = 3  # number of groups

# Between-group sum of squares (SSB)
ssb = (len(group_a) * (group_a.mean() - grand_mean)**2 +
       len(group_b) * (group_b.mean() - grand_mean)**2 +
       len(group_c) * (group_c.mean() - grand_mean)**2)

# Within-group sum of squares (SSW)
ssw = (np.sum((group_a - group_a.mean())**2) +
       np.sum((group_b - group_b.mean())**2) +
       np.sum((group_c - group_c.mean())**2))

# Degrees of freedom
df_between = k - 1
df_within = n_total - k

# Mean squares
msb = ssb / df_between
msw = ssw / df_within

# F-statistic
f_statistic_manual = msb / msw

# P-value from F distribution
p_value_manual = stats.f.sf(f_statistic_manual, df_between, df_within)

print("Manual ANOVA Calculation:")
print(f"  SSB: {ssb:.4f}, SSW: {ssw:.4f}")
print(f"  MSB: {msb:.4f}, MSW: {msw:.4f}")
print(f"  F-statistic: {f_statistic_manual:.4f}")
print(f"  P-value: {p_value_manual:.6f}")

# Verify with scipy.stats.f_oneway
f_statistic_scipy, p_value_scipy = stats.f_oneway(group_a, group_b, group_c)
print(f"\nSciPy f_oneway:")
print(f"  F-statistic: {f_statistic_scipy:.4f}")
print(f"  P-value: {p_value_scipy:.6f}")

# Interpretation
alpha = 0.05
critical_value = stats.f.ppf(1 - alpha, df_between, df_within)
print(f"\nCritical value (α=0.05): {critical_value:.4f}")

if p_value_scipy < alpha:
    print("Result: Significant difference between group means")
else:
    print("Result: No significant difference between group means")

The manual calculation matches SciPy’s result, demonstrating how the F distribution underlies ANOVA. In production code, use f_oneway for convenience and reliability.

Visualization and Summary Statistics

Understanding how degrees of freedom affect the F distribution’s shape helps you interpret results correctly.

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

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Panel 1: Effect of numerator df
x = np.linspace(0.01, 5, 500)
dfd_fixed = 30
for dfn in [1, 3, 10, 30]:
    axes[0, 0].plot(x, stats.f.pdf(x, dfn, dfd_fixed), 
                    label=f'dfn={dfn}', linewidth=2)
axes[0, 0].set_title(f'Varying Numerator df (dfd={dfd_fixed})')
axes[0, 0].set_xlabel('x')
axes[0, 0].set_ylabel('Density')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# Panel 2: Effect of denominator df
dfn_fixed = 5
for dfd in [5, 10, 30, 100]:
    axes[0, 1].plot(x, stats.f.pdf(x, dfn_fixed, dfd), 
                    label=f'dfd={dfd}', linewidth=2)
axes[0, 1].set_title(f'Varying Denominator df (dfn={dfn_fixed})')
axes[0, 1].set_xlabel('x')
axes[0, 1].set_ylabel('Density')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# Panel 3: Summary statistics table
df_values = [(2, 10), (5, 20), (10, 30), (20, 50), (50, 100)]
stats_data = []

for dfn, dfd in df_values:
    f_dist = stats.f(dfn, dfd)
    mean = f_dist.mean() if dfd > 2 else np.nan
    var = f_dist.var() if dfd > 4 else np.nan
    skew = f_dist.stats(moments='s')[0] if dfd > 6 else np.nan
    stats_data.append([dfn, dfd, mean, var, skew])

axes[1, 0].axis('off')
table = axes[1, 0].table(
    cellText=[[f"{v:.3f}" if isinstance(v, float) else str(v) for v in row] 
              for row in stats_data],
    colLabels=['dfn', 'dfd', 'Mean', 'Variance', 'Skewness'],
    loc='center',
    cellLoc='center'
)
table.auto_set_font_size(False)
table.set_fontsize(10)
table.scale(1.2, 1.5)
axes[1, 0].set_title('Summary Statistics for Various df Combinations')

# Panel 4: Critical values at different significance levels
alphas = [0.10, 0.05, 0.01]
dfn_range = range(1, 31)
dfd_fixed = 20

for alpha in alphas:
    critical_vals = [stats.f.ppf(1 - alpha, dfn, dfd_fixed) for dfn in dfn_range]
    axes[1, 1].plot(dfn_range, critical_vals, label=f'α={alpha}', linewidth=2)

axes[1, 1].set_title(f'Critical Values vs Numerator df (dfd={dfd_fixed})')
axes[1, 1].set_xlabel('Numerator df')
axes[1, 1].set_ylabel('Critical Value')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Notice how increasing degrees of freedom shifts the distribution toward symmetry and reduces critical values. This directly impacts statistical power—larger samples give you more precise variance estimates and lower thresholds for detecting real effects.

The F distribution is a workhorse of statistical inference. Master these Python tools, understand the underlying mathematics, and you’ll handle variance comparisons and ANOVA with confidence.