29 - Quanto aspettare il terzo successo? La distribuzione binomiale negativa

Nell'articolo 14 abbiamo lanciato un dado finché non è comparso il primo 6.

La domanda era:

Quanti lanci servono per ottenere il primo successo?

La distribuzione geometrica risponde esattamente a questo problema.

Ora cambiamo appena la domanda:

Quanti lanci servono per ottenere il terzo 6?

Non fissiamo più in anticipo il numero di lanci. Continuiamo finché raggiungiamo un numero prefissato di successi.

Questa è la situazione descritta dalla distribuzione binomiale negativa.

Binomiale e binomiale negativa guardano lo stesso processo da due lati

Con una sequenza di prove indipendenti, tutte con la stessa probabilità di successo p, possiamo porre due domande diverse.

La binomiale chiede:

fisso n prove
→ quanti successi ottengo?

La binomiale negativa chiede:

fisso r successi
→ quante prove servono?

Nel nostro esempio:

successo = esce 6
p = 1/6
r = 3

Chiamiamo:

X = numero totale di lanci fino al terzo 6

Il valore minimo possibile è:

$$ X = 3 $$

perché possiamo ottenere:

6 6 6

al primo, secondo e terzo lancio.

Non esiste invece un massimo finito. Potremmo attendere molto a lungo prima di osservare il terzo 6.

Quindi il supporto della distribuzione — cioè l'insieme dei valori che X può assumere con probabilità positiva — è:

$$ X = 3, 4, 5, ... $$

Una domanda concreta: il terzo 6 arriva esattamente al decimo lancio

Vogliamo calcolare:

$$ P(X=10) $$

Se il terzo 6 arriva esattamente al decimo lancio, il decimo lancio deve essere necessariamente:

Nei primi nove lanci devono invece esserci esattamente due 6.

Le loro posizioni possono essere scelte in:

$$ \binom{9}{2} = 36 $$

modi.

Per una particolare disposizione dei primi nove lanci abbiamo:

2 successi
7 insuccessi

quindi la probabilità è:

$$ (1/6)^2 (5/6)^7 $$

Il decimo lancio deve poi essere un altro successo:

$$ 1/6 $$

Perciò:

$$ P(X=10) = \binom{9}{2} (1/6)^3 (5/6)^7 $$

Il risultato esatto è:

$$ 78125 / 1679616 $$

cioè circa:

4,65136%

Perché compare C(n-1,r-1)

Generalizziamo.

Vogliamo che il successo numero r avvenga esattamente alla prova n.

La prova n deve quindi essere un successo.

Prima di essa, nelle prime:

$$ n-1 $$

prove, dobbiamo avere esattamente:

$$ r-1 $$

successi.

Le loro posizioni possono essere scelte in:

$$ \binom{n-1}{r-1} $$

modi.

In totale avremo:

r successi
n-r insuccessi

La probabilità è quindi:

$$ P(X=n) = \binom{n-1}{r-1} p^r (1-p)^{n-r} $$

per:

$$ n = r, r+1, r+2, ... $$

Questa è la funzione di massa della binomiale negativa nella convenzione che useremo nella serie.

Attenzione: esistono due convenzioni

Il nome binomiale negativa viene usato con due variabili leggermente diverse.

Una convenzione conta:

Y = numero di insuccessi prima del successo numero r

quindi:

$$ Y = 0,1,2,... $$

La nostra convenzione conta invece:

X = numero totale di prove fino al successo numero r

Le due variabili sono semplicemente collegate da:

$$ X = Y + r $$

Per evitare ambiguità, in tutta questa serie useremo X come numero totale di prove.

È anche la convenzione più naturale rispetto all'articolo 14, dove il tempo geometrico comprendeva il lancio del primo successo.

La geometrica è il caso r=1

Mettiamo:

$$ r = 1 $$

nella formula:

$$ P(X=n) = \binom{n-1}{0} p (1-p)^{n-1} $$

ma:

$$ \binom{n-1}{0}=1 $$

quindi:

$$ P(X=n) = p(1-p)^{n-1} $$

che è esattamente la distribuzione geometrica già incontrata.

La binomiale negativa non sostituisce quindi la geometrica: la generalizza.

geometrica
→ attesa fino al primo successo

binomiale negativa
→ attesa fino al successo numero r

Un altro collegamento: la stessa domanda vista con la binomiale

Supponiamo di chiedere:

Qual è la probabilità che il terzo 6 sia già arrivato entro 10 lanci?

In termini della nostra variabile:

$$ P(X \le 10) $$

Ma dire che il terzo 6 è già arrivato entro il decimo lancio equivale a dire:

nei primi 10 lanci sono comparsi almeno tre 6.

Se:

B ~ Binomiale(10,1/6)

allora:

$$ P(X \le 10) = P(B \ge 3) $$

Calcolando la coda binomiale otteniamo:

P(X <= 10)
=
566299 / 2519424
≈
22,4773%

Questo collegamento è importante:

binomiale negativa
→ tempo necessario per raggiungere r successi

binomiale
→ successi accumulati entro n prove

Le due descrizioni parlano dello stesso processo di Bernoulli.

Quanto dobbiamo aspettare in media?

Possiamo interpretare il tempo fino al terzo successo come somma di tre attese.

Chiamiamo:

G1 = prove fino al primo successo
G2 = prove dal primo al secondo successo
G3 = prove dal secondo al terzo successo

Allora:

$$ X = G1 + G2 + G3 $$

Ogni Gi ha distribuzione geometrica con la stessa probabilità p.

Per una geometrica che conta le prove fino al successo:

$$ E[G] = 1/p $$

Quindi, per linearità del valore atteso:

$$ E[X] = r/p $$

Nel nostro caso:

E[X]
=
3 / (1/6)
=
18

Servono dunque 18 lanci in media per osservare il terzo 6.

Questo non significa che il terzo 6 arriverà normalmente proprio al lancio 18.

La distribuzione è piuttosto dispersa e ha una lunga coda verso destra.

Per esempio:

P(X <= 18)
≈
59,7346%

Quindi, anche arrivati al numero medio di lanci, resta una probabilità superiore al 40% che il terzo 6 non sia ancora comparso.

Anche la varianza nasce dalla somma delle attese

Per una geometrica nella nostra convenzione:

Var(G)
=
(1-p)/p^2

Le attese fra successi successivi sono indipendenti nel modello di prove di Bernoulli indipendenti.

Possiamo quindi sommare le varianze:

$$ \operatorname{Var}(X) = r(1-p)/p^2 $$

Per il terzo 6:

Var(X)
=
3 · (5/6) / (1/6)^2
=
90

La deviazione standard è:

sqrt(90)
≈
9,4868 lanci

La media 18 va quindi letta insieme alla dispersione: il tempo di attesa può oscillare parecchio da una sequenza all'altra.

Un caso più piccolo: terza testa con una moneta equa

Se il successo è testa:

$$ \begin{gathered} p = 1/2 \ r = 3 \end{gathered} $$

allora:

E[X]
=
3/(1/2)
=
6

ed:

Var(X)
=
3(1/2)/(1/2)^2
=
6

Le prime probabilità sono:

$$ \begin{gathered} P(X=3) = 1/8 \ P(X=4) = 3/16 \ P(X=5) = 3/16 \ P(X=6) = 5/32 \end{gathered} $$

Sono già sufficienti per vedere che la distribuzione non assegna tutta la massa attorno alla media.

Il calcolo in C#

Lo standalone C# associato a questo articolo è:

BinomialeNegativa.cs

Possiamo implementare direttamente la funzione di massa.

using System.Numerics;

static BigInteger Combinazioni(int n, int k)
{
    if (n < 0)
        throw new ArgumentOutOfRangeException(nameof(n));

    if (k < 0 || k > n)
        return BigInteger.Zero;

    k = Math.Min(k, n - k);

    BigInteger risultato = BigInteger.One;

    for (int i = 1; i <= k; i++)
    {
        risultato *= n - k + i;
        risultato /= i;
    }

    return risultato;
}

static double BinomialeNegativa(
    int proveTotali,
    int successiRichiesti,
    double probabilitaSuccesso)
{
    if (successiRichiesti <= 0)
        throw new ArgumentOutOfRangeException(
            nameof(successiRichiesti));

    if (probabilitaSuccesso <= 0.0 ||
        probabilitaSuccesso > 1.0)
    {
        throw new ArgumentOutOfRangeException(
            nameof(probabilitaSuccesso));
    }

    if (proveTotali < successiRichiesti)
        return 0.0;

    return
        (double)Combinazioni(
            proveTotali - 1,
            successiRichiesti - 1)
        *
        Math.Pow(
            probabilitaSuccesso,
            successiRichiesti)
        *
        Math.Pow(
            1.0 - probabilitaSuccesso,
            proveTotali - successiRichiesti);
}

Per il terzo 6 esattamente al decimo lancio:

double p = BinomialeNegativa(
    proveTotali: 10,
    successiRichiesti: 3,
    probabilitaSuccesso: 1.0 / 6.0);

Console.WriteLine($"{p:P9}");

otteniamo circa:

4,651360787%

Il file standalone associato all'articolo aggiunge:

confronto con la geometrica quando r=1;
coda cumulativa;
verifica della relazione con la binomiale;
media e varianza teoriche;
simulazione a seed fisso come controllo secondario.

Le ipotesi del modello contano

La binomiale negativa richiede prove:

indipendenti
con due esiti rispetto all'evento studiato
con probabilità di successo p costante

Per il dado ideale, ogni lancio mantiene:

$$ p(6)=1/6 $$

anche dopo cento lanci senza 6.

Se invece estraiamo carte senza reinserimento, la probabilità di successo cambia dopo ogni carta osservata.

In quel caso il modello binomiale negativo non è più corretto.

Questa distinzione sarà precisamente il punto di partenza del prossimo articolo.

Il punto pratico

La binomiale e la binomiale negativa non descrivono due fenomeni separati.

Descrivono due domande sullo stesso processo:

Binomiale:
quanti successi ottengo in n prove?

Binomiale negativa:
quante prove servono per ottenere r successi?

Per la nostra convenzione:

$$ P(X=n) = \binom{n-1}{r-1} p^r (1-p)^{n-r} $$

con:

$$ \begin{gathered} E[X] = r/p \ \operatorname{Var}(X) = r(1-p)/p^2 \end{gathered} $$

La geometrica ricompare come caso particolare r=1.

Nel prossimo articolo manterremo la domanda "quanto devo aspettare il successo numero r?", ma toglieremo l'ipotesi di reinserimento: lavoreremo con una popolazione finita in cui la composizione cambia a ogni estrazione. Nascerà così l'ipergeometrica inversa.