spb/traite-evaluation-qc

1"""Chapitre 39 — Modèles bayésiens et intervalles : vérification des identités, des diagnostics et des chiffres cités dans le texte,2les tableaux et les solutions. Tolérances déclarées dans chaque assertion. Tout est numpy/scipy (pymc et arviz absents)."""3import sys, pathlib, math4import numpy as np, pytest5sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1] / "cases"))6sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))7import ch39_bayes as c398from traite import cas1910@pytest.fixture(scope="module")11def res():12    return c39.main(verbose=False)1314# ----------------------------------------------------------------------------------------------------------------- conjugué15def test_nig_reference_coincide_avec_mco(res):16    ex = res["conjugue"]; p = ex["p600"]; m = res["modele"]17    assert np.allclose(p["m_n"], m["beta"], atol=1e-9)                      # a priori de référence : m_n = β̂18    assert p["a_n"] == pytest.approx(292.0) and p["nu"] == pytest.approx(584.0)   # (n − k)/2 = (600 − 16)/219    assert p["b_n"] * 2 == pytest.approx(float(m["residus"] @ m["residus"]), rel=1e-6)   # b_n = SCR/2 (différence d'arrondi y'y − m'V⁻¹m)20    assert p["sigma_moy"] == pytest.approx(0.0832, abs=0.0002)21    t = ex["tab600"]; assert int(t["couvre"].sum()) == 15 and int(t["vrai"].notna().sum()) == 1522    assert t.loc["ln_S", "post_sd"] == pytest.approx(0.0138, abs=0.0002)23    assert t.loc["sect_Masson-Angers", "post_moy"] == pytest.approx(-0.2631, abs=0.0005)2425def test_prior_propre_vague(res):26    ex = res["conjugue"]; pp = ex["p600_propre"]27    assert np.abs(pp["m_n"] - res["modele"]["beta"]).max() < 1e-428    assert pp["a_n"] == pytest.approx(300.001, abs=1e-6) and pp["sigma_moy"] == pytest.approx(0.0821, abs=0.0002)29    assert ex["pred600_propre"]["bas95"] == pytest.approx(540_678, abs=5) and ex["pred600_propre"]["haut95"] == pytest.approx(747_707, abs=5)3031def test_predictive_student_vs_ic_frequentiste(res):32    pr = res["conjugue"]["pred600"]33    assert pr["mediane"] == pytest.approx(635_821, abs=2)34    assert pr["echelle"] == pytest.approx(0.08364, abs=0.00005) and pr["part_coefficients"] * 100 == pytest.approx(1.49, abs=0.02)35    assert pr["bas95"] == pytest.approx(539_504, abs=5) and pr["haut95"] == pytest.approx(749_334, abs=5)36    assert pr["bas90"] == pytest.approx(553_979, abs=5) and pr["haut90"] == pytest.approx(729_754, abs=5)37    # écart avec l'intervalle fréquentiste du ch. 18 (HC1, quantile normal) inférieur à 0,1 %38    assert abs(pr["bas95"] / c39.IC18["bas95"] - 1) < 0.001 and abs(pr["haut95"] / c39.IC18["haut95"] - 1) < 0.00139    # identité : quantile retransformé = exp(loc ± t·échelle)40    from scipy import stats41    q = stats.t.ppf(0.975, pr["nu"]) * pr["echelle"]42    assert math.exp(pr["loc"] - q) == pytest.approx(pr["bas95"], rel=1e-12)4344def test_coherence_sequentielle(res):45    ex = res["conjugue"]; p = ex["p600"]; ps = ex["p_rec_seq"]46    assert ex["n_anc"] + ex["n_rec"] == 60047    assert np.abs(ps["m_n"] - p["m_n"]).max() < 1e-8          # postérieure(2023-24) comme a priori + ventes 2025-26 = postérieure(600)48    assert ps["a_n"] == pytest.approx(p["a_n"], abs=1e-9) and ps["b_n"] == pytest.approx(p["b_n"], rel=1e-6)49    assert np.allclose(ps["V_n"], p["V_n"], atol=1e-10)5051def test_petit_echantillon_et_prior_informatif(res):52    ex = res["conjugue"]; cols = ex["cols_micro"]; i = cols.index("ln_S")53    assert ex["n_micro"] == 40 and len(cols) == 1254    f = ex["petits"]["faible"]; k1 = ex["petits"]["kappa_1"]; k4 = ex["petits"]["kappa_4"]; k16 = ex["petits"]["kappa_16"]55    assert f["post"]["m_n"][i] == pytest.approx(0.611, abs=0.002) and f["post"]["sd"][i] == pytest.approx(0.067, abs=0.002)56    assert k1["post"]["m_n"][i] == pytest.approx(0.543, abs=0.002) and k1["post"]["sd"][i] == pytest.approx(0.016, abs=0.001)57    assert k4["post"]["sd"][i] == pytest.approx(0.028, abs=0.001) and k16["post"]["sd"][i] == pytest.approx(0.042, abs=0.001)58    assert f["pred"]["mediane"] == pytest.approx(650_965, abs=5) and k1["pred"]["mediane"] == pytest.approx(630_978, abs=5)59    assert f["largeur95_pct"] == pytest.approx(39.5, abs=0.1) and k1["largeur95_pct"] == pytest.approx(33.2, abs=0.1)60    # plus κ est grand, moins l'a priori pèse : écart-type croissant, moyenne qui se rapproche de l'estimation plate61    assert k1["post"]["sd"][i] < k4["post"]["sd"][i] < k16["post"]["sd"][i] < f["post"]["sd"][i]62    assert ex["pred_hull"]["mediane"] == pytest.approx(623_894, abs=5) and ex["p_hull"]["n"] == 1866364def test_loo_exact_nig_petit():65    """LOO exact : sur un petit jeu simulé, la prédictive LOO calculée par réestimation égale celle d'une réestimation explicite."""66    rng = np.random.default_rng(1); X = np.column_stack([np.ones(30), rng.normal(size=30)]); y = 1 + 0.5 * X[:, 1] + rng.normal(0, 0.3, 30)67    loo = c39.loo_exact_nig(X, y, c39.prior_reference(["const", "x"]))68    from scipy import stats69    p = c39.nig_posterieure(X[1:], y[1:], c39.prior_reference(["const", "x"]))70    s2 = p["b_n"] / p["a_n"] * (1 + X[0] @ p["V_n"] @ X[0])71    assert loo["pointwise"][0] == pytest.approx(stats.t.logpdf(y[0], p["nu"], loc=X[0] @ p["m_n"], scale=math.sqrt(s2)), rel=1e-10)7273# ----------------------------------------------------------------------------------------------------------------- hiérarchique74def test_gibbs_diagnostics(res):75    d = res["hier"]["diag"]76    assert (d["rhat_bloc"] < 1.01).all() and (d["ess_bloc"] > 2000).all()77    assert d.loc["dist_centre", "ess_naif"] < 300 and d.loc["dist_centre", "rhat_naif"] > 1.03   # colinéarité distance/secteurs : mélange lent78    assert d.loc["alpha_Masson-Angers", "ess_naif"] < 300 and d.loc["alpha_Masson-Angers", "ess_bloc"] > 40007980def test_hierarchique_retrouve_les_lois(res):81    hi = res["hier"]; r = hi["resume"]; rel = hi["relatifs"]82    assert r.loc["sigma", "moyenne"] == pytest.approx(0.0832, abs=0.0015)83    assert r.loc["ln_S", "moyenne"] == pytest.approx(0.5395, abs=0.003)84    for s in rel.index:   # chaque prime de secteur vraie est dans l'intervalle de crédibilité à 95 %85        assert rel.loc[s, "q025"] <= rel.loc[s, "vrai"] <= rel.loc[s, "q975"]86    assert rel.loc["Masson-Angers", "moyenne"] == pytest.approx(-0.224, abs=0.01)87    assert math.sqrt(hi["tau2_med"]) == pytest.approx(0.18, abs=0.02)88    # rétrécissement : l'effet aléatoire est moins extrême que l'effet fixe pour les secteurs périphériques89    for s in ["Buckingham", "Masson-Angers"]:90        assert rel.loc[s, "moyenne"] > rel.loc[s, "fixe"]9192def test_formule_retrecissement_et_experience(res):93    hi = res["hier"]; e = hi["experience"]94    assert c39.facteur_retrecissement(5, 0.083 ** 2, 0.15 ** 2) == pytest.approx(0.9423, abs=0.0005)95    assert c39.facteur_retrecissement(5, 0.083 ** 2, 0.05 ** 2) == pytest.approx(0.6447, abs=0.0005)96    assert list(e.index) == [5, 10, 20, 47]97    a0_rel = float(hi["g"]["alpha0"].mean() - hi["g"]["alpha"][:, :, 0].mean())98    for n in e.index:99        assert e.loc[n, "hier"] > e.loc[n, "fixe"]                                   # rétrécissement vers la moyenne des secteurs100        assert abs(e.loc[n, "hier"] - a0_rel) < abs(e.loc[n, "fixe"] - a0_rel)101        assert abs(e.loc[n, "hier"] - e.loc[n, "vrai"]) < abs(e.loc[n, "fixe"] - e.loc[n, "vrai"])   # ici, plus proche du vrai102        assert e.loc[n, "hier_rhat"] < 1.01103    assert e.loc[5, "fixe"] == pytest.approx(-0.342, abs=0.003) and e.loc[5, "hier"] == pytest.approx(-0.300, abs=0.01)104    assert e.loc[5, "B"] == pytest.approx(0.970, abs=0.005)105    assert e.loc[5, "largeur_f_pct"] == pytest.approx(36.1, abs=0.3) and e.loc[47, "largeur_f_pct"] == pytest.approx(33.4, abs=0.3)106    # la part du bruit de transaction domine : la largeur ne varie que de quelques points entre 5 et 47 ventes107    assert e.loc[5, "largeur_h_pct"] - e.loc[47, "largeur_h_pct"] < 5108109def test_predictive_hierarchique(res):110    pr = res["hier"]["pred"]111    assert pr["mediane"] == pytest.approx(634_748, rel=0.004) and pr["moyenne"] == pytest.approx(637_727, rel=0.004)112    assert pr["bas95"] == pytest.approx(540_309, rel=0.006) and pr["haut95"] == pytest.approx(749_431, rel=0.006)113    assert pr["sd_log"] == pytest.approx(0.0834, abs=0.002) and pr["sd_moyenne_log"] == pytest.approx(0.0103, abs=0.002)114    # la moyenne a posteriori en dollars est proche de la prédiction de Duan du ch. 18 (637 953 $) : la retransformation est implicite115    assert abs(pr["moyenne"] / c39.IC18["point"] - 1) < 0.005116117def test_sensibilite_prior_tau(res):118    s = res["hier"]["sensibilite"]119    assert len(s) == 4 and s["tau_med_600"].between(0.12, 0.20).all()120    assert (s["ma_5"].max() - s["ma_5"].min()) < 0.05          # l'effet à 5 ventes bouge de moins de 0,05 selon l'a priori sur τ121    assert (s["ma_600"].max() - s["ma_600"].min()) < 0.03122123def test_criteres(res):124    c = res["criteres"]125    pool, fixe, hier = c.index126    assert c.loc[fixe, "elpd_waic"] - c.loc[pool, "elpd_waic"] > 100127    assert abs(c.loc[fixe, "elpd_waic"] - c.loc[hier, "elpd_waic"]) < 2128    for m in (pool, fixe):   # WAIC ≈ LOO exact ≈ LOO tronqué pour les modèles conjugués129        assert abs(c.loc[m, "elpd_waic"] - c.loc[m, "elpd_loo_exact"]) < 1.0130        assert abs(c.loc[m, "elpd_loo_tronque"] - c.loc[m, "elpd_loo_exact"]) < 0.5131    assert 10 < c.loc[fixe, "p_waic"] < 20 and c.loc[fixe, "waic"] == pytest.approx(-1266.3, abs=1.0)132133# ----------------------------------------------------------------------------------------------------------------- simulation134def test_recuperation_sur_donnees_simulees(res):135    """Données engendrées par les lois du générateur (README synthétique) : le modèle conjugué couvre ≥ 13 des 15 paramètres et le136    modèle hiérarchique retrouve σ et les primes de secteur dans ses intervalles de crédibilité."""137    v = res["ventes"]; d = c39.simuler_depuis_generateur(v, 600, graine=3)138    p = c39.ajuster_nig(d); cols = p["colonnes"]139    vrais = np.array([c39.VRAIS.get(c, np.nan) for c in cols]); m = ~np.isnan(vrais)140    couv = ((vrais[m] >= p["ic95_bas"][m]) & (vrais[m] <= p["ic95_haut"][m])).sum()141    assert couv >= 13 and p["sigma_moy"] == pytest.approx(0.085, abs=0.007)142    Xc, moy, z = c39.matrice_hier(d); y = np.log(d["prix_vente"].astype(float)).values143    g = c39.gibbs_hierarchique(y, Xc.values, z, n_iter=1500, brulage=500, chaines=2, graine=11)144    rel = c39.effets_relatifs(g); r = c39.resume_gibbs(g, list(Xc.columns))145    assert r.loc["sigma", "moyenne"] == pytest.approx(0.085, abs=0.007)146    ok = sum(rel.loc[s, "q025"] <= c39.VRAIS[f"sect_{s}"] <= rel.loc[s, "q975"] for s in rel.index)147    assert ok >= 3 and (r["rhat"] < 1.05).all()148149def test_rhat_et_ess_sur_chaines_independantes():150    rng = np.random.default_rng(3); ch = rng.normal(size=(4, 2000))151    assert 0.99 < c39.rhat_scinde(ch) < 1.02 and c39.ess(ch) > 0.7 * 8000152    ar = np.empty((4, 2000)); ar[:, 0] = rng.normal(size=4)      # AR(1) à φ = 0,9 : ESS ≈ mn (1−φ)/(1+φ) ≈ 420153    for t in range(1, 2000): ar[:, t] = 0.9 * ar[:, t - 1] + rng.normal(size=4) * math.sqrt(1 - 0.81)154    assert 200 < c39.ess(ar) < 800 and c39.rhat_scinde(ar) < 1.1155    assert c39.rhat_scinde(ar + np.array([[0.0], [0.0], [3.0], [0.0]])) > 1.5   # une chaîne décalée : R̂ explose156157def test_waic_identites():158    rng = np.random.default_rng(5); lp = rng.normal(-1, 0.1, size=(500, 50))159    w = c39.waic(lp)160    assert w["elpd_waic"] == pytest.approx(w["lppd"] - w["p_waic"]) and w["waic"] == pytest.approx(-2 * w["elpd_waic"])161    assert w["p_waic"] == pytest.approx(50 * 0.01, rel=0.2)162163# ----------------------------------------------------------------------------------------------------------------- exercices164def test_exercices(res):165    e = res["exercices"]166    assert e["ex1"]["m_n"] == pytest.approx(3440.8, abs=0.1) and e["ex1"]["s_n"] == pytest.approx(120.0, abs=0.01)167    assert e["ex1"]["poids_donnees"] == pytest.approx(0.64, abs=0.001) and e["ex1"]["n_equivalent_prior"] == pytest.approx(2.25)168    assert e["ex2"]["p_sup_650"] == pytest.approx(0.396, abs=0.002) and e["ex2"]["p_inf_550"] == pytest.approx(0.042, abs=0.002)169    assert e["ex3"][5] == pytest.approx(0.9423, abs=0.0005) and e["ex3"]["alpha_5"] == pytest.approx(-0.2873, abs=0.0005)170    assert e["ex3"]["sd_cond_5"] == pytest.approx(0.0360, abs=0.0002) and e["ex3"]["sd_fixe_5"] == pytest.approx(0.0371, abs=0.0002)171    assert e["ex4"]["t975"] == pytest.approx(1.9640, abs=0.0002) and e["ex4"]["t975_20"] == pytest.approx(2.0860, abs=0.0002)172    assert e["ex5"]["rhat"] == pytest.approx(1.0077, abs=0.0002) and e["ex5"]["rhat_divergent"] == pytest.approx(1.2855, abs=0.001)173    assert e["ex7"]["m13.0_s1.0"]["sd"] == pytest.approx(0.081, abs=0.003) and e["ex7"]["m13.0_s0.05"]["sd"] == pytest.approx(0.0425, abs=0.002)174