"""04. Potencia (militar: los controles de 2021 no cuentan, no hay positivas del compuesto ese año) (Hanley-McNeil 1982) con los recuentos de 03 (positivas y controles del compuesto). Python puro."""
import math
N = {  # cártel: (positivas del compuesto, controles del compuesto, licitaciones positivas, órganos positivos)
    'alimentos S/0016/21': (87, 1318, 33, 15), 'militar S/0008/21': (20, 708, 15, 10), 'carreteras S/0013/19': (4, 10, 4, 1),
    'conjunto (3)': (111, 2036, 52, 26)}
def se_hm(A, n1, n0):
    q1 = A / (2 - A); q2 = 2 * A * A / (1 + A)
    return math.sqrt((A * (1 - A) + (n1 - 1) * (q1 - A * A) + (n0 - 1) * (q2 - A * A)) / (n1 * n0))
def phi(z): return 0.5 * (1 + math.erf(z / math.sqrt(2)))
out = ["| Cártel | positivas / controles (licitaciones; órganos) | deff | semiancho IC95 con A=0,65 | P(AUC^ ≥ 0,65) si A = 0,50 / 0,60 / 0,65 / 0,70 / 0,80 | P(LB95 > 0,50) si A = 0,65 / 0,70 / 0,80 |",
       "|---|---|---|---|---|---|"]
for k, (n1, n0, nl, no) in N.items():
    for deff in (1, 2, 4):
        s65 = se_hm(0.65, n1, n0) * math.sqrt(deff)
        pp = [1 - phi((0.65 - A) / (se_hm(A, n1, n0) * math.sqrt(deff))) for A in (0.5, 0.6, 0.65, 0.7, 0.8)]
        pl = [1 - phi((0.5 + 1.96 * se_hm(A, n1, n0) * math.sqrt(deff) - A) / (se_hm(A, n1, n0) * math.sqrt(deff))) for A in (0.65, 0.7, 0.8)]
        out.append(f"| {k} | {n1} / {n0} ({nl}; {no}) | {deff} | ±{1.96*s65:.3f} | " + " / ".join(f"{x:.2f}" for x in pp) + " | " + " / ".join(f"{x:.2f}" for x in pl) + " |")
# probabilidad de «2 de 3» por azar (A = 0,5 en los tres), deff 2
p = [1 - phi((0.65 - 0.5) / (se_hm(0.5, n1, n0) * math.sqrt(2))) for (n1, n0, _, _) in list(N.values())[:3]]
p2 = p[0]*p[1]*(1-p[2]) + p[0]*(1-p[1])*p[2] + (1-p[0])*p[1]*p[2] + p[0]*p[1]*p[2]
out.append(f"\nFalso «funciona» por azar con un solo candidato (A = 0,50 en los tres, deff 2): P(≥ 2 de 3 con AUC^ ≥ 0,65) = {p2:.3f} "
           f"(por cártel {', '.join(f'{x:.3f}' for x in p)}). Carreteras pasa el umbral por azar con probabilidad {p[2]:.2f}.")
print("\n".join(out))