"""Chapitre 39 — Modèles bayésiens et intervalles : vérification des identités, des diagnostics et des chiffres cités dans le texte, les tableaux et les solutions. Tolérances déclarées dans chaque assertion. Tout est numpy/scipy (pymc et arviz absents).""" import sys, pathlib, math import numpy as np, pytest sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1] / "cases")) sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1])) import ch39_bayes as c39 from traite import cas1 @pytest.fixture(scope="module") def res(): return c39.main(verbose=False) # ----------------------------------------------------------------------------------------------------------------- conjugué def test_nig_reference_coincide_avec_mco(res): ex = res["conjugue"]; p = ex["p600"]; m = res["modele"] assert np.allclose(p["m_n"], m["beta"], atol=1e-9) # a priori de référence : m_n = β̂ assert p["a_n"] == pytest.approx(292.0) and p["nu"] == pytest.approx(584.0) # (n − k)/2 = (600 − 16)/2 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) assert p["sigma_moy"] == pytest.approx(0.0832, abs=0.0002) t = ex["tab600"]; assert int(t["couvre"].sum()) == 15 and int(t["vrai"].notna().sum()) == 15 assert t.loc["ln_S", "post_sd"] == pytest.approx(0.0138, abs=0.0002) assert t.loc["sect_Masson-Angers", "post_moy"] == pytest.approx(-0.2631, abs=0.0005) def test_prior_propre_vague(res): ex = res["conjugue"]; pp = ex["p600_propre"] assert np.abs(pp["m_n"] - res["modele"]["beta"]).max() < 1e-4 assert pp["a_n"] == pytest.approx(300.001, abs=1e-6) and pp["sigma_moy"] == pytest.approx(0.0821, abs=0.0002) assert ex["pred600_propre"]["bas95"] == pytest.approx(540_678, abs=5) and ex["pred600_propre"]["haut95"] == pytest.approx(747_707, abs=5) def test_predictive_student_vs_ic_frequentiste(res): pr = res["conjugue"]["pred600"] assert pr["mediane"] == pytest.approx(635_821, abs=2) assert pr["echelle"] == pytest.approx(0.08364, abs=0.00005) and pr["part_coefficients"] * 100 == pytest.approx(1.49, abs=0.02) assert pr["bas95"] == pytest.approx(539_504, abs=5) and pr["haut95"] == pytest.approx(749_334, abs=5) assert pr["bas90"] == pytest.approx(553_979, abs=5) and pr["haut90"] == pytest.approx(729_754, abs=5) # écart avec l'intervalle fréquentiste du ch. 18 (HC1, quantile normal) inférieur à 0,1 % assert abs(pr["bas95"] / c39.IC18["bas95"] - 1) < 0.001 and abs(pr["haut95"] / c39.IC18["haut95"] - 1) < 0.001 # identité : quantile retransformé = exp(loc ± t·échelle) from scipy import stats q = stats.t.ppf(0.975, pr["nu"]) * pr["echelle"] assert math.exp(pr["loc"] - q) == pytest.approx(pr["bas95"], rel=1e-12) def test_coherence_sequentielle(res): ex = res["conjugue"]; p = ex["p600"]; ps = ex["p_rec_seq"] assert ex["n_anc"] + ex["n_rec"] == 600 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) assert ps["a_n"] == pytest.approx(p["a_n"], abs=1e-9) and ps["b_n"] == pytest.approx(p["b_n"], rel=1e-6) assert np.allclose(ps["V_n"], p["V_n"], atol=1e-10) def test_petit_echantillon_et_prior_informatif(res): ex = res["conjugue"]; cols = ex["cols_micro"]; i = cols.index("ln_S") assert ex["n_micro"] == 40 and len(cols) == 12 f = ex["petits"]["faible"]; k1 = ex["petits"]["kappa_1"]; k4 = ex["petits"]["kappa_4"]; k16 = ex["petits"]["kappa_16"] 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) 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) assert k4["post"]["sd"][i] == pytest.approx(0.028, abs=0.001) and k16["post"]["sd"][i] == pytest.approx(0.042, abs=0.001) assert f["pred"]["mediane"] == pytest.approx(650_965, abs=5) and k1["pred"]["mediane"] == pytest.approx(630_978, abs=5) assert f["largeur95_pct"] == pytest.approx(39.5, abs=0.1) and k1["largeur95_pct"] == pytest.approx(33.2, abs=0.1) # plus κ est grand, moins l'a priori pèse : écart-type croissant, moyenne qui se rapproche de l'estimation plate assert k1["post"]["sd"][i] < k4["post"]["sd"][i] < k16["post"]["sd"][i] < f["post"]["sd"][i] assert ex["pred_hull"]["mediane"] == pytest.approx(623_894, abs=5) and ex["p_hull"]["n"] == 186 def test_loo_exact_nig_petit(): """LOO exact : sur un petit jeu simulé, la prédictive LOO calculée par réestimation égale celle d'une réestimation explicite.""" 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) loo = c39.loo_exact_nig(X, y, c39.prior_reference(["const", "x"])) from scipy import stats p = c39.nig_posterieure(X[1:], y[1:], c39.prior_reference(["const", "x"])) s2 = p["b_n"] / p["a_n"] * (1 + X[0] @ p["V_n"] @ X[0]) 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) # ----------------------------------------------------------------------------------------------------------------- hiérarchique def test_gibbs_diagnostics(res): d = res["hier"]["diag"] assert (d["rhat_bloc"] < 1.01).all() and (d["ess_bloc"] > 2000).all() assert d.loc["dist_centre", "ess_naif"] < 300 and d.loc["dist_centre", "rhat_naif"] > 1.03 # colinéarité distance/secteurs : mélange lent assert d.loc["alpha_Masson-Angers", "ess_naif"] < 300 and d.loc["alpha_Masson-Angers", "ess_bloc"] > 4000 def test_hierarchique_retrouve_les_lois(res): hi = res["hier"]; r = hi["resume"]; rel = hi["relatifs"] assert r.loc["sigma", "moyenne"] == pytest.approx(0.0832, abs=0.0015) assert r.loc["ln_S", "moyenne"] == pytest.approx(0.5395, abs=0.003) for s in rel.index: # chaque prime de secteur vraie est dans l'intervalle de crédibilité à 95 % assert rel.loc[s, "q025"] <= rel.loc[s, "vrai"] <= rel.loc[s, "q975"] assert rel.loc["Masson-Angers", "moyenne"] == pytest.approx(-0.224, abs=0.01) assert math.sqrt(hi["tau2_med"]) == pytest.approx(0.18, abs=0.02) # rétrécissement : l'effet aléatoire est moins extrême que l'effet fixe pour les secteurs périphériques for s in ["Buckingham", "Masson-Angers"]: assert rel.loc[s, "moyenne"] > rel.loc[s, "fixe"] def test_formule_retrecissement_et_experience(res): hi = res["hier"]; e = hi["experience"] assert c39.facteur_retrecissement(5, 0.083 ** 2, 0.15 ** 2) == pytest.approx(0.9423, abs=0.0005) assert c39.facteur_retrecissement(5, 0.083 ** 2, 0.05 ** 2) == pytest.approx(0.6447, abs=0.0005) assert list(e.index) == [5, 10, 20, 47] a0_rel = float(hi["g"]["alpha0"].mean() - hi["g"]["alpha"][:, :, 0].mean()) for n in e.index: assert e.loc[n, "hier"] > e.loc[n, "fixe"] # rétrécissement vers la moyenne des secteurs assert abs(e.loc[n, "hier"] - a0_rel) < abs(e.loc[n, "fixe"] - a0_rel) assert abs(e.loc[n, "hier"] - e.loc[n, "vrai"]) < abs(e.loc[n, "fixe"] - e.loc[n, "vrai"]) # ici, plus proche du vrai assert e.loc[n, "hier_rhat"] < 1.01 assert e.loc[5, "fixe"] == pytest.approx(-0.342, abs=0.003) and e.loc[5, "hier"] == pytest.approx(-0.300, abs=0.01) assert e.loc[5, "B"] == pytest.approx(0.970, abs=0.005) 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) # la part du bruit de transaction domine : la largeur ne varie que de quelques points entre 5 et 47 ventes assert e.loc[5, "largeur_h_pct"] - e.loc[47, "largeur_h_pct"] < 5 def test_predictive_hierarchique(res): pr = res["hier"]["pred"] assert pr["mediane"] == pytest.approx(634_748, rel=0.004) and pr["moyenne"] == pytest.approx(637_727, rel=0.004) assert pr["bas95"] == pytest.approx(540_309, rel=0.006) and pr["haut95"] == pytest.approx(749_431, rel=0.006) assert pr["sd_log"] == pytest.approx(0.0834, abs=0.002) and pr["sd_moyenne_log"] == pytest.approx(0.0103, abs=0.002) # la moyenne a posteriori en dollars est proche de la prédiction de Duan du ch. 18 (637 953 $) : la retransformation est implicite assert abs(pr["moyenne"] / c39.IC18["point"] - 1) < 0.005 def test_sensibilite_prior_tau(res): s = res["hier"]["sensibilite"] assert len(s) == 4 and s["tau_med_600"].between(0.12, 0.20).all() 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 τ assert (s["ma_600"].max() - s["ma_600"].min()) < 0.03 def test_criteres(res): c = res["criteres"] pool, fixe, hier = c.index assert c.loc[fixe, "elpd_waic"] - c.loc[pool, "elpd_waic"] > 100 assert abs(c.loc[fixe, "elpd_waic"] - c.loc[hier, "elpd_waic"]) < 2 for m in (pool, fixe): # WAIC ≈ LOO exact ≈ LOO tronqué pour les modèles conjugués assert abs(c.loc[m, "elpd_waic"] - c.loc[m, "elpd_loo_exact"]) < 1.0 assert abs(c.loc[m, "elpd_loo_tronque"] - c.loc[m, "elpd_loo_exact"]) < 0.5 assert 10 < c.loc[fixe, "p_waic"] < 20 and c.loc[fixe, "waic"] == pytest.approx(-1266.3, abs=1.0) # ----------------------------------------------------------------------------------------------------------------- simulation def test_recuperation_sur_donnees_simulees(res): """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 le modèle hiérarchique retrouve σ et les primes de secteur dans ses intervalles de crédibilité.""" v = res["ventes"]; d = c39.simuler_depuis_generateur(v, 600, graine=3) p = c39.ajuster_nig(d); cols = p["colonnes"] vrais = np.array([c39.VRAIS.get(c, np.nan) for c in cols]); m = ~np.isnan(vrais) couv = ((vrais[m] >= p["ic95_bas"][m]) & (vrais[m] <= p["ic95_haut"][m])).sum() assert couv >= 13 and p["sigma_moy"] == pytest.approx(0.085, abs=0.007) Xc, moy, z = c39.matrice_hier(d); y = np.log(d["prix_vente"].astype(float)).values g = c39.gibbs_hierarchique(y, Xc.values, z, n_iter=1500, brulage=500, chaines=2, graine=11) rel = c39.effets_relatifs(g); r = c39.resume_gibbs(g, list(Xc.columns)) assert r.loc["sigma", "moyenne"] == pytest.approx(0.085, abs=0.007) ok = sum(rel.loc[s, "q025"] <= c39.VRAIS[f"sect_{s}"] <= rel.loc[s, "q975"] for s in rel.index) assert ok >= 3 and (r["rhat"] < 1.05).all() def test_rhat_et_ess_sur_chaines_independantes(): rng = np.random.default_rng(3); ch = rng.normal(size=(4, 2000)) assert 0.99 < c39.rhat_scinde(ch) < 1.02 and c39.ess(ch) > 0.7 * 8000 ar = np.empty((4, 2000)); ar[:, 0] = rng.normal(size=4) # AR(1) à φ = 0,9 : ESS ≈ mn (1−φ)/(1+φ) ≈ 420 for t in range(1, 2000): ar[:, t] = 0.9 * ar[:, t - 1] + rng.normal(size=4) * math.sqrt(1 - 0.81) assert 200 < c39.ess(ar) < 800 and c39.rhat_scinde(ar) < 1.1 assert c39.rhat_scinde(ar + np.array([[0.0], [0.0], [3.0], [0.0]])) > 1.5 # une chaîne décalée : R̂ explose def test_waic_identites(): rng = np.random.default_rng(5); lp = rng.normal(-1, 0.1, size=(500, 50)) w = c39.waic(lp) assert w["elpd_waic"] == pytest.approx(w["lppd"] - w["p_waic"]) and w["waic"] == pytest.approx(-2 * w["elpd_waic"]) assert w["p_waic"] == pytest.approx(50 * 0.01, rel=0.2) # ----------------------------------------------------------------------------------------------------------------- exercices def test_exercices(res): e = res["exercices"] assert e["ex1"]["m_n"] == pytest.approx(3440.8, abs=0.1) and e["ex1"]["s_n"] == pytest.approx(120.0, abs=0.01) assert e["ex1"]["poids_donnees"] == pytest.approx(0.64, abs=0.001) and e["ex1"]["n_equivalent_prior"] == pytest.approx(2.25) 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) assert e["ex3"][5] == pytest.approx(0.9423, abs=0.0005) and e["ex3"]["alpha_5"] == pytest.approx(-0.2873, abs=0.0005) 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) assert e["ex4"]["t975"] == pytest.approx(1.9640, abs=0.0002) and e["ex4"]["t975_20"] == pytest.approx(2.0860, abs=0.0002) assert e["ex5"]["rhat"] == pytest.approx(1.0077, abs=0.0002) and e["ex5"]["rhat_divergent"] == pytest.approx(1.2855, abs=0.001) 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)