Cheatsheet : scipy.stats

scipy.stats est le module Python de référence pour les distributions de probabilité, les tests statistiques et l'estimation.

Distributions de probabilité

Créer et manipuler une loi

from scipy import stats

# Loi normale μ=2, σ²=4
X = stats.norm(loc=2, scale=2)  # scale = écart-type σ

# Densité de probabilité (PDF)
X.pdf(3)           # f(3)

# Fonction de répartition (CDF)
X.cdf(2)           # P(X ≤ 2)

# Quantile (inverse de la CDF)
X.ppf(0.975)       # x tel que P(X ≤ x) = 0.975

# Espérance et variance
X.mean(), X.var(), X.std()

# Générer un échantillon
echantillon = X.rvs(size=1000, random_state=42)

Lois continues courantes

# Normale
X = stats.norm(loc=mu, scale=sigma)

# Uniforme sur [a, b]
X = stats.uniform(loc=a, scale=b - a)

# Exponentielle de paramètre λ
X = stats.expon(scale=1/lam)

# Student à ν degrés de liberté
X = stats.t(df=nu)

# Chi-deux à ν degrés de liberté
X = stats.chi2(df=nu)

Lois discrètes courantes

# Bernoulli
X = stats.bernoulli(p=0.3)

# Binomiale
X = stats.binom(n=20, p=0.5)

# Poisson
X = stats.poisson(mu=5)

Loi normale multivariée

import numpy as np

mu = np.array([1, 2])
Sigma = np.array([[2, 0.5],
                   [0.5, 1]])

X = stats.multivariate_normal(mean=mu, cov=Sigma)

X.pdf([1, 2])                      # densité en un point
echantillon = X.rvs(size=500)      # échantillon 2D

Statistiques descriptives

Sur un échantillon

import numpy as np

x = stats.norm(loc=5, scale=2).rvs(size=100, random_state=42)

# Statistiques de base
np.mean(x)          # moyenne empirique
np.var(x, ddof=1)   # variance corrigée (ddof=1)
np.std(x, ddof=1)   # écart-type corrigé
np.median(x)        # médiane

# Statistiques descriptives complètes
stats.describe(x)

Moments

# Moments centrés d'ordre k
stats.moment(x, moment=3)   # moment d'ordre 3

# Asymétrie (skewness) et aplatissement (kurtosis)
stats.skew(x)
stats.kurtosis(x)

Estimation

Maximum de vraisemblance (fit)

# Ajuster une loi normale aux données
mu_hat, sigma_hat = stats.norm.fit(x)
print(f"μ̂ = {mu_hat:.3f}, σ̂ = {sigma_hat:.3f}")

# Ajuster une loi exponentielle
loc_hat, scale_hat = stats.expon.fit(x)
lam_hat = 1 / scale_hat

# Ajuster en fixant certains paramètres
mu_hat, sigma_hat = stats.norm.fit(x, floc=0)  # fixer μ=0

Intervalles de confiance

# IC pour la moyenne (variance inconnue, loi de Student)
n = len(x)
x_bar = np.mean(x)
s = np.std(x, ddof=1)
alpha = 0.05

# Méthode manuelle
t_crit = stats.t.ppf(1 - alpha/2, df=n-1)
ic = (x_bar - t_crit * s / np.sqrt(n),
      x_bar + t_crit * s / np.sqrt(n))

# Méthode directe
ic = stats.t.interval(1 - alpha, df=n-1,
                       loc=x_bar, scale=s/np.sqrt(n))
print(f"IC à 95% : [{ic[0]:.3f}, {ic[1]:.3f}]")

Borne de Cramér-Rao (exemple gaussien)

# Information de Fisher pour N(μ, σ²) par rapport à μ
# I(μ) = 1/σ², donc CRB = σ²/n
sigma2 = 4
n = 50
crb = sigma2 / n
print(f"Borne de Cramér-Rao : {crb:.4f}")

# Vérification : variance de la moyenne empirique
N_sim = 10000
estimations = [np.mean(stats.norm(0, np.sqrt(sigma2)).rvs(n))
               for _ in range(N_sim)]
print(f"Variance empirique de μ̂ : {np.var(estimations):.4f}")

Tests statistiques

Test de normalité

# Test de Shapiro-Wilk
stat, p_value = stats.shapiro(x)
print(f"Shapiro-Wilk : stat={stat:.4f}, p={p_value:.4f}")

# Test de Kolmogorov-Smirnov
stat, p_value = stats.kstest(x, 'norm', args=(np.mean(x), np.std(x)))
print(f"KS : stat={stat:.4f}, p={p_value:.4f}")

Test de Student (comparaison de moyennes)

# Test sur un échantillon : H₀: μ = μ₀
stat, p_value = stats.ttest_1samp(x, popmean=5)

# Test sur deux échantillons indépendants : H₀: μ₁ = μ₂
x1 = stats.norm(5, 1).rvs(50, random_state=1)
x2 = stats.norm(5.5, 1).rvs(50, random_state=2)
stat, p_value = stats.ttest_ind(x1, x2)

# Test sur échantillons appariés
stat, p_value = stats.ttest_rel(x1, x2)

Test du Chi-deux

# Test d'adéquation
observed = np.array([18, 22, 20, 15, 25])
expected = np.array([20, 20, 20, 20, 20])
stat, p_value = stats.chisquare(observed, f_exp=expected)

Recettes courantes

Simuler un estimateur et tracer sa distribution

import numpy as np
from scipy import stats

# Paramètres
mu_vrai, sigma_vrai = 3, 2
n = 30
N_sim = 5000

# Simulations
estimations = np.array([
    np.mean(stats.norm(mu_vrai, sigma_vrai).rvs(n))
    for _ in range(N_sim)
])

# Vérification des propriétés
print(f"E[μ̂] = {np.mean(estimations):.3f}  (attendu : {mu_vrai})")
print(f"Var(μ̂) = {np.var(estimations):.4f}  (attendu : {sigma_vrai**2/n:.4f})")

Comparaison loi théorique vs empirique (QQ-plot)

from scipy import stats

# QQ-plot
fig, ax = plt.subplots()
stats.probplot(x, dist="norm", plot=ax)
ax.set_title("QQ-plot")

Log-vraisemblance

# Calcul de la log-vraisemblance pour une loi normale
def log_vraisemblance(mu, sigma, x):
    return np.sum(stats.norm.logpdf(x, loc=mu, scale=sigma))

# Grille de μ pour tracer la log-vraisemblance
mu_grid = np.linspace(0, 10, 200)
ll = [log_vraisemblance(mu, sigma_vrai, x) for mu in mu_grid]