L11 — Quante osservazioni servono per distinguere un vantaggio dal caso

from math import ceil, erf, sqrt

import matplotlib.pyplot as plt
import numpy as np

from cvbook import seed_for
from cvbook.dati import carica
from cvbook.metriche import rendimenti

VANTAGGIO = 0.001        # ← guadagno medio per operazione, al netto dei costi
                         # PROVA / TRY: il TUO vantaggio stimato (esercizio 1)
OSCILLAZIONE = 0.035     # ← deviazione standard del risultato per operazione
                         # PROVA / TRY: la TUA oscillazione (esercizio 1)
OPERAZIONI_ANNO = 250    # ← quante ne fai in un anno
POTENZA = 0.80           # ← probabilita' di accorgersene, se il vantaggio esiste
ALFA = 0.05              # ← rischio accettato di scambiare rumore per segnale


def quantile_normale(p: float) -> float:
    """Inversa della normale standard, per bisezione: nessuna dipendenza esterna."""
    basso, alto = -10.0, 10.0
    for _ in range(200):
        mezzo = (basso + alto) / 2
        if 0.5 * (1 + erf(mezzo / sqrt(2))) < p:
            basso = mezzo
        else:
            alto = mezzo
    return (basso + alto) / 2


def quante_servono(vantaggio: float, oscillazione: float,
                   potenza: float = POTENZA, alfa: float = ALFA) -> int:
    if vantaggio <= 0:
        raise ValueError("il vantaggio deve essere positivo")
    z_alfa = quantile_normale(1 - alfa)
    z_potenza = quantile_normale(potenza)
    return int(ceil(((z_alfa + z_potenza) * oscillazione / vantaggio) ** 2))


n = quante_servono(VANTAGGIO, OSCILLAZIONE)
print(f"vantaggio da dimostrare: {VANTAGGIO:.3%} per operazione")
print(f"oscillazione per operazione: {OSCILLAZIONE:.1%}\n")
print(f"operazioni necessarie: {n:,}")
print(f"a {OPERAZIONI_ANNO} operazioni l'anno: {n / OPERAZIONI_ANNO:,.1f} anni")

oscillazione_btc = float(np.std(rendimenti(
    carica("btcusdt").sort("data")["chiusura"].to_numpy()), ddof=1))
print(f"oscillazione giornaliera misurata su Bitcoin: {oscillazione_btc:.1%}\n")

vantaggi = [0.0005, 0.001, 0.002, 0.005, 0.010]
print(f"{'vantaggio':>10s} {'operazioni':>12s} {'a 250/anno':>14s}")
for v in vantaggi:
    q = quante_servono(v, oscillazione_btc)
    print(f"{v:10.2%} {q:12,d} {q / 250:13.1f} anni")

with avvio.figura("schermo"):
    fig, ax = plt.subplots()
    griglia = np.linspace(0.0003, 0.012, 200)
    ax.plot(griglia * 100, [quante_servono(v, oscillazione_btc) for v in griglia], linewidth=2)
    for v in vantaggi:
        ax.plot([v * 100], [quante_servono(v, oscillazione_btc)], marker="o")
    ax.set_yscale("log")
    ax.set_xlabel("Vantaggio medio per operazione (%)")
    ax.set_ylabel("Operazioni necessarie (scala log)")
    plt.show()

print("\nUn vantaggio dello 0,1% per operazione sarebbe un risultato eccellente — "
      "e i costi si mangiano gia' lo 0,12% a giro. Per dimostrarlo servono "
      "decenni di operativita' quotidiana.")

LANCI = 40
VITTORIE = 26
PROVE = 200_000  # PROVA / TRY: 20000 (veloce) · 200000 (percentuale più precisa)

rng = np.random.default_rng(seed_for("lab-potere-moneta"))
esiti = rng.binomial(LANCI, 0.5, PROVE)
quota = float((esiti >= VITTORIE).mean())

print(f"su {PROVE:,} sequenze di {LANCI} lanci di una moneta EQUA:")
print(f"  {VITTORIE} vittorie o piu': {quota:.2%} delle volte")
print(f"\nE se hai provato piu' di una manciata di strategie prima di trovare "
      f"questa, quel {quota:.0%} te lo sei praticamente garantito.")

STRATEGIE = 20        # PROVA / TRY: 200 (vedi esercizio 2)
OSSERVAZIONI = 500    # PROVA / TRY: 5000 (vedi esercizio 3)

rng = np.random.default_rng(seed_for("lab-potere-multipli"))
soglia = quantile_normale(1 - ALFA)

print(f"{'strategia':>10s} {'risultato medio':>17s} {'statistica t':>14s} {'supera il test?':>17s}")
passate = 0
for k in range(STRATEGIE):
    campione = rng.normal(0.0, OSCILLAZIONE, OSSERVAZIONI)  # vantaggio ESATTAMENTE zero
    t = campione.mean() / (campione.std(ddof=1) / np.sqrt(OSSERVAZIONI))
    supera = t > soglia
    passate += supera
    print(f"{k + 1:10d} {campione.mean():17.4%} {t:14.2f} {'SI' if supera else 'no':>17s}")

print(f"\n{passate} strategie su {STRATEGIE} hanno superato il test standard.")
print("Dentro non c'era niente. Nessuna di esse aveva un vantaggio: era zero, "
      "messo li' da noi.")

def soglia_corretta(tentativi: int, alfa: float = ALFA) -> float:
    """Correzione conservativa: si divide il rischio accettato per i tentativi."""
    return quantile_normale(1 - alfa / tentativi)


print(f"{'tentativi':>10s} {'soglia sulla statistica t':>27s} "
      f"{'prob. che almeno uno passi':>28s}")
for tentativi in (1, 5, 20, 50, 100, 500):
    print(f"{tentativi:10d} {soglia_corretta(tentativi):27.2f} "
          f"{1 - (1 - ALFA) ** tentativi:27.1%}")

print("\nCon cento tentativi, trovare qualcosa che supera il test non corretto e' "
      "praticamente certo. Non e' un difetto del test: e' la sua definizione.")