"""Monte Carlo for the Messi-Yamal chain posterior.

Model (all conditioning at the 2007 information state I; later data used
only for time-stable rates with protagonists excluded -- the Evidence Rule).
p1..p5 below are the five factors of the post's chain rule:
p1=Pr(E_win), p2=Pr(E_pair|E_win), p3=Pr(E_Yamal), p4=Pr(E_Messi),
p5=Pr(E_final|E_Yamal, E_Messi), with E_win the raffle win, E_pair the
pairing with Messi, E_Yamal Yamal on Spain's 2026 squad, E_Messi Messi
on Argentina's, E_final an Argentina-Spain final.

  p1: raffle win. N entrants ~ LogNormal, elicited: 90% mass in [30, 1000].
      One winner drawn. No data; prior = posterior. p1 = 1/N.
  p2: paired with Messi. Prior K ~ Uniform{10..14} over calendar
      player-pages, updated on the documented count of 12 monthly pages
      -> point mass at K=12. Conflicting published rosters are about who
      filled the pages, not how many there were; Messi is taken to occupy
      exactly one. Exchangeable assignment => p2 = 1/12 exactly.
  p3: boy born 2007 reaches Spain's 2026 squad aged 18-19.
      lambda = teenage (<=19) squad members per Spain tournament PLAYED,
      so Spain's 2026 qualification is treated as certain (8 straight
      appearances as of 2007; pricing it at 0.9 moves the mean ~10%).
      Prior Gamma(1, 2) (structural, mean 0.5).
      Data: k=5 teenagers in n=17 tournaments (Yamal excluded).
      Posterior Gamma(6, 19). The <=19 window spans ~2 birth cohorts of
      ~254,000 boys each. 2007-knowable family/geography boost
      m ~ LogNormal(ln 2, 0.35). p3 = lambda * m / (2 * 254_000).
  p4: established 20-year-old outfield star still in a WC squad 20 years on.
      Rate marginalizes over the nation qualifying: a cohort star whose
      nation missed 2026 counts as a non-survivor. Prior Beta(0.5, 0.5)
      (Jeffreys).
      Data: 2006 U21 outfield cohort, Messi excluded: k=1 of n=23 (Modric).
      Posterior Beta(1.5, 22.5). p4 = q.
  p5: Argentina and Spain meet in the final, given both elite.
      r = P(a pre-tournament top-favorite side reaches a given final).
      Prior Beta(1, 2). Data: pre-tournament favorite reached the final
      4 of 15 times, 1966-2022 (SportsOddsHistory/covers.com series).
      Structure factor 1/2 for opposite halves. p5 = 0.5 * r_A * r_B
      (two independent draws of r, teams are exchangeable draws from the
      elite class).
"""
import numpy as np

rng = np.random.default_rng(20260719)
S = 2_000_000

# p1
mu_N, sd_N = np.log(np.sqrt(30 * 1000)), np.log(1000 / np.sqrt(30 * 1000)) / 1.645
N = rng.lognormal(mu_N, sd_N, S)
p1 = 1.0 / N

# p2
p2 = np.full(S, 1.0 / 12.0)

# p3
lam = rng.gamma(1 + 5, 1 / (2 + 17), S)
m = rng.lognormal(np.log(2), 0.35, S)
COHORT_BOYS = 254_000
p3 = lam * m / (2 * COHORT_BOYS)

# p4
q = rng.beta(0.5 + 1, 0.5 + 22, S)
p4 = q

# p5
k5, n5 = 4, 15
rA = rng.beta(1 + k5, 2 + n5 - k5, S)
rB = rng.beta(1 + k5, 2 + n5 - k5, S)
p5 = 0.5 * rA * rB

P = p1 * p2 * p3 * p4 * p5

def one_in(x):
    return f"1 in {1/x:,.0f}"

print("factor means:")
for name, p in [("p1", p1), ("p2", p2), ("p3", p3), ("p4", p4), ("p5", p5)]:
    print(f"  {name}: mean={p.mean():.3e}  ({one_in(p.mean())})")
print()
print(f"posterior mean of P : {P.mean():.3e}  ({one_in(P.mean())})")
print(f"posterior median    : {np.median(P):.3e}  ({one_in(np.median(P))})")
q5, q95 = np.quantile(P, [0.05, 0.95])
print(f"90% credible interval: [{q5:.3e}, {q95:.3e}]")
print(f"  i.e. between {one_in(q95)} and {one_in(q5)}")