Problema n.8

Testo del problema

Un individuo può avere un certo numero di discendenti diretti: con probabilità (Subscript[p, 1])ha un discendente, con probabilità (Subscript[p, 2])ha due discendenti, con probabilità (Subscript[p, 0])non ha alcun discendente ((Subscript[p, 0]) + (Subscript[p, 1]) + (Subscript[p, 2]) = 1).
Ciascuno dei discendenti a sua volta dà luogo a propri discendenti (seconda generazione) con la stessa legge di probabilità, e così di seguito nelle successive generazioni.
Studiare la probabilità di estinzione del processo. (Un modello matematico di questo genere si applica allo studio dei neutroni nei processi a cascata).

La soluzione esatta

Indichiamo con P(n) la probabilità che n individui si estinguano "prima o poi". Possiamo immaginare di suddividere l'insieme di n individui in due sottoinsiemi disgiunti, che evidentemente non si influenzano a vicenda per quanto riguarda la probabilità di estinzione; l'insieme complessivo di individui si estingue se si estinguono entrambi i gruppi in cui è stato diviso. Si ha perciò la regola:

         P(m + n)  =  P(m) P(n)

ed in particolare:

         P(n + 1)  =  P(n) P(1),

da cui si ricava subito:

         P(n)  =  ((P(1))^n).

A noi però interessa proprio P(1); per determinarlo, osserviamo che un singolo individuo può dare origine, nella generazione successiva, a 0, 1 o 2 individui, con probabilità rispettive (Subscript[p, 0]), (Subscript[p, 1]) e (Subscript[p, 2]), quindi possiamo scrivere:

         \(TraditionalForm
      \`P(1)\  = p\_0\  + \ p\_1\ \(P(1)\)\  + p\_2\ \(P(2)\)\),

cioè:

        \(TraditionalForm
      \`P(1)\  = p\_0\  + \ p\_1\ \(P(1)\)\  + p\_2\ \(P(1)\)\^2\).

Questa è un'equazione di secondo grado in P(1) (se (Subscript[p, 2]) > 0) e tenendo conto che (Subscript[p, 0]) + (Subscript[p, 1]) + (Subscript[p, 2]) = 1, troviamo le due soluzioni:

         P(1)  = (Subscript[p, 0])/(Subscript[p, 2])    e    P(1)  =  1.

Per  (Subscript[p, 0]) -> 0 la probabilità di estinzione si deve annullare e quindi  per (Subscript[p, 0]) vicino a zero la soluzione corretta è la prima: P(1)  = (Subscript[p, 0])/(Subscript[p, 2]). Tale soluzione non ha però senso per (Subscript[p, 0]) > (Subscript[p, 2]) quindi in tale caso bisogna prendere la seconda soluzione:  P(1)  =  1. 

I casi particolari si trattano agevolmente: se (Subscript[p, 2]) = 0  si ha  P(1)  =  1 finché (Subscript[p, 0]) > 0, altrimenti  P(1)  =  0, come ovvio, se (Subscript[p, 0]) = (Subscript[p, 2]) = 0.

Approssimazioni successive

La soluzione trovata nella sezione precedente non richiede l'uso di un software particolare per essere trovata. Possiamo però chiederci, più in generale, quale sia la probabilità p(k, n) di  avere  n individui dopo k generazioni; la probabilità di estinzione dopo k generazioni sarà quindi p(k, 0)  e per k -> ∞ dovremmo ritrovare la probabilità "asintotica" trovata in precedenza.

Vogliamo trovare una procedura iterativa per calcolare p(k + 1, n) quando sia nota p(k, n) per ogni n. Cominciamo con l'osservare che all'inizio (generazione zero) si ha: p(0, 1) = 1 e p(0, n)  =  0 se n ≠ 1. Inoltre alla generazione k-sima si hanno al massimo (2^k) individui.

Supponiamo ora di conoscere tutte le probabilità p(k, n)  della k-sima generazione, per n che va da 0 a (2^k). Se abbiamo n individui, di questi (Subscript[n, 0])avranno 0 discendenti, (Subscript[n, 1])avranno 1 discendente e (Subscript[n, 2]) avranno 2 discendenti (con (Subscript[n, 0]) + (Subscript[n, 1]) + (Subscript[n, 2]) = n), quindi avremo (Subscript[n, 1]) + 2  (Subscript[n, 2]) individui nella generazione successiva, e questo con probabilità 

       (({{n}, {(Subscript[n, 0])}}))(Subscript[(Subscript[p, 0]^n), 0])(({{n - (Subscript[n, 0])}, {(Subscript[n, 1])}}))(Subscript[(Subscript[p, 1]^n), 1])(Subscript[(Subscript[p, 2]^n), 2]). 

Non è difficile capire il motivo per cui sono presenti i coefficienti binomiali: gli (Subscript[n, 0]) individui che non avranno discendenti si possono scegliere in (({{n}, {(Subscript[n, 0])}})) modi, col che rimangono n - (Subscript[n, 0])indivudui fra cui scegliere gli (Subscript[n, 1]) che avranno un discendente, ecc. Questa probabilità, moltiplicata naturalmente per p(k, n), va ad incrementare dunque la probabilità p(k + 1, (Subscript[n, 1]) + 2  (Subscript[n, 2])) di avere (Subscript[n, 1]) + 2  (Subscript[n, 2]) individui nella generazione successiva; questo vale per tutti i possibili modi di suddividere n in (Subscript[n, 0]) + (Subscript[n, 1]) + (Subscript[n, 2]) e per n che va da 0 a (2^k).

Proviamo ad implementare questo metodo con Mathematica. Cominciamo col dare le definizioni preliminari e poniamo uguali a zero tutti i p(k, n) eccetto p(0, 1) = 1.

\(Clear[p]; 
p[_, \ _]\  := \ 0; 
p[0, \ 1]\  = \ 1; \)

La funzione compgen consiste di un ciclo Do quadruplo: k numera le successive generazioni, n indica il numero di individui alla generazione k-sima e (Subscript[n, 0]), (Subscript[n, 1]), (Subscript[n, 2]) quanti fra gli n individui hanno avuto 0, 1 o 2 discendenti. (Abbiamo però usato il fatto che (Subscript[n, 0]) + (Subscript[n, 1]) + (Subscript[n, 2]) = n per eliminare (Subscript[n, 2]).) Alla fine del ciclo tutti i valori di p(k, n) saranno stati calcolati fino a k = ngen.

\(compgen[ngen_, p0_, p1_, p2_] := 
		
      Do[p[k + 1, 2 \((n - n0)\) - n1]\  += \ 
		\ \ \ \ \ \ \ \ \ 
          p[k, \ n]*\(p0\^n0\) Binomial[n, n0]*

            			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            p1\^n1\ Binomial[n - n0, n1]*

            			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            p2\^\(n - n0 - n1\), 

        	\ \ \ \ \ \ {k, 0, ngen - 1}, {n, \ 0, 2\^k\ }, {n0, 0, n}, {n1, \ 
          0, \ n - n0}]\)

Il tempo di elaborazione cresce però molto rapidamente all'aumentare di ngen, sicché ci limitiamo alla quinta generazione. Notare che le probabilità (Subscript[p, 0]), (Subscript[p, 1]) e (Subscript[p, 2]) sono lasciate in forma simbolica, in modo che il risultato sarà un polinomio in queste tre variabili.

compgen[5, p0, p1, p2]

Ecco i risultati ottenuti per le probabilità di estinzione; come si vede, p(k, n) è un polinomio di grado (2^k) - 1 in (Subscript[p, 0]), (Subscript[p, 1]) e (Subscript[p, 2]).

p[1, 0] // Expand

p[2, 0] // Expand

p[3, 0] // Expand

p[4, 0] // Expand

p[5, 0] // Expand

Vogliamo ora visualizzare i risultati ottenuti e confrontarli con la soluzione asintotica trovata prima. Il metodo migliore sembra essere quello di fissare un valore per (Subscript[p, 1])e fare quindi un grafico di p(k, 0) al variare di (Subscript[p, 0]) fra 0 e 1 - (Subscript[p, 1]). Allo stesso grafico sovrapponiamo (in rosso) il grafico della funzione limite plim:

\(plim := \ p0/p2\  /; \ p0 < p2; 

    plim\  := \ 1\  /; \ p0 ≥ p2; \)

Cominciamo col caso p1=0.

p1 = 0;  p2 =  1 - p0;

Costruiamo i grafici per i valori di k da 1 a 5, usando un ciclo Do. I 5 grafici risultanti possono essere visti in successione; come si vede, p(k, 0) converge abbastanza rapidamente al valore limite.

\(Do[

      Plot[{plim, p[k, 0]}, \ {p0, 0, 1 - p1}, 
	
        PlotStyle -> {Hue[1], GrayLevel[0]}, 
	PlotRange -> All], 
      
{k, 1, 5}]\)

Se facciamo la stessa cosa cambiando il valore di p1 non emergono differenze significative, salvo il fatto che la convergenza è un po' più lenta.

p1 = 0.3;  p2 =  1 - p0 - p1;

\(Do[

      Plot[{plim, p[k, 0]}, \ {p0, 0, 1 - p1}, 
	
        PlotStyle -> {Hue[1], GrayLevel[0]}, 
	PlotRange -> All], 
      
{k, 1, 5}]\)

Se ci limitiamo al caso particolare (Subscript[p, 1]) = 0 allora possiamo semplificare notevolmente il calcolo di cui sopra ed arrivare agevolmente fino alla settima generazione.

\(Clear[p, p0, p1, p2]; 
p[_, \ _]\  := \ 0; 
p[0, \ 1]\  = \ 1; 

    k\  = \ 0; \)

\(\(	Do[
      p[k + 1, 2  n - 2  n0]\  += \ 
		\ \ \ \ \ \ \ \ \ 
        p[k, \ n]*\(p0\^n0\) Binomial[n, n0] \((1 - p0)\)\^\(n - n0\), 

      	\ \ \ \ \ \ {n, \ 0, 2\^k\ }, {n0, 0, n}]; 
\(++k\)\)\)

Facciamo ora gli stessi grafici di prima; come si vede, all'aumentare di k il grafico di p(k, 0) si avvicina sempre più al risultato asintotico tracciato in rosso.

\(Do[

      Plot[{plim, p[k, 0]}, \ {p0, 0, 1}, 
	
        PlotStyle -> {Hue[1], GrayLevel[0]}, 
	PlotRange -> All], 
      
{k, 1, 7}]\)

Uso del metodo Monte-Carlo

Ci proponiamo ora di verificare la soluzione asintotica mediante una simulazione che usa la generazione di numeri pseudo-casuali, ottenuti in Mathematica tramite la funzione Random.

La funzione figli[(Subscript[p, 0]),  (Subscript[p, 1])] trova per il singolo individuo in modo casuale il numero di discendenti (0, 1 o 2) con probabilità (Subscript[p, 0]),  (Subscript[p, 1]) e (Subscript[p, 2]) = 1 - (Subscript[p, 0]) - (Subscript[p, 1]). Questo viene fatto usando la funzione Random[], che genera un numero pseudo-casuale compreso fra 0 e 1 con distribuzione uniforme; questo numero appartiene quindi all'intervallo [0, (Subscript[p, 0])] con probabilità (Subscript[p, 0]), mentre è nell'intervallo [(Subscript[p, 0]), (Subscript[p, 0]) + (Subscript[p, 1])] con probabilità (Subscript[p, 1]) ed infine appartiene a  [(Subscript[p, 0]) + (Subscript[p, 1]), 1] con probabilità(Subscript[p, 2]). La scelta fra queste alternative viene effettuata con la funzione Which.

Clear["Global`*"]

\(figli[p0_, p1_] := 
	
      Module[{r = Random[]}, 
		
        Which[r <= p0, \ 0, \ r\  <= \ \((p0 + \ p1)\), \ 1, \ True, \ 2]\ \ ]
        \)

Come esempio, simuliamo la discendenza avuta da 100 individui nella generazione successiva, quando (Subscript[p, 0]) = (Subscript[p, 1]) = (Subscript[p, 2]) = 1/3.

Table[figli[1/3, 1/3],  {100}]

Se sommiamo tutti i numeri di questa lista, otteniamo il numero di discendenti nella generazione successiva:

Plus @@ %

Possiamo fare la stessa cosa usando la funzione Sum. Notare che il risultato è diverso da prima, perché diversi saranno, in generale, i numeri casuali generati nei due casi.

Sum[figli[1/3, 1/3],  {100}]

Iterando questa procedura possiamo generare in modo casuale una successione di generazioni. Nell'esempio che segue, partiamo con 20 individui e li seguiamo per 10 generazioni.

NestList[Sum[figli[1/3, 1/3], {#}]&, 20,  10]

Definiamo dunque un'altra funzione, nuovagen[p0,p1,n], che a partire da una popolazione di n individui trova in modo casuale il numero di individui nella generazione successiva. È conveniente definirla come funzione compilata, per renderne più veloce l'esecuzione.

nuovagen=\(Compile[{p0, p1, {n, _Integer}}, 
		\ 
        Module[{r, p = p0 + p1}, 
			
          Sum[	Which[\((r = Random[])\) < p0, 0, r\  < p, 1, True, 2], 

            				{n}]
		\ \ ]]\)

Ecco come appare lo stesso esempio di prima, usando questa nuova funzione.

NestList[nuovagen[1/3, 1/3, #]&, 20,  10]

Possiamo confrontare le velocità di elaborazione nei due casi. Abbiamo usato SeedRandom per inizializzare la generazione dei numeri casuali nello stesso modo ed ottenere dunque gli stessi risultati.

\(SeedRandom[1]; 

    Timing[NestList[Sum[figli[1./3, 1./3], {#}]&, 20, \ 20]]\)

\(SeedRandom[1]; 

    Timing[NestList[nuovagen[1./3, 1./3, #]&, 20, \ 20]]\)

Per valutare empiricamente la probabilità di estinzione, possiamo procedere così: partendo con 1 individuo iteriamo nuovagen, finché otteniamo zero oppure finché il numero di individui supera un certo valore limite nsoglia al di sopra della quale possiamo ritenere la probabilità di estinzione trascurabile. Ripetiamo questa procedura un numero niter (sufficientemente elevato) di volte; il numero di volte che abbiamo ottenuto zero, diviso per il numero totale di prove, ci darà una stima della probabilità. Questa stima sarà tanto migliore quanto maggiori sono i valori di nsoglia e di niter; il valore ottimale va determinato per tentativi e risulterà da un compromesso fra l'accuratezza dei risultati e la velocità dell'esecuzione.

Proviamo per esempio con questi valori:

\(p0 = \ 0.2; \ p2 = 0.4; \ p1 = 1 - p0 - p2; \ 
nsoglia\  = \ 10; \ 
    niter = 100; \)

Ecco come appare la procedura: il ciclo While continua ad iterare nuovagen (partendo con n=1) finché otteniamo zero oppure un numero di individui maggiore di nsoglia; il ciclo Do ripete il tutto un numero niter di volte. La variabile nzeri viene incrementata di una unità ogni volta che nuovagen dà come risultato zero, quindi conta il numero di estinzioni che si sono verificate.

\(\(n = 1; \ nzeri = 0; 
\ 
    Do[\ While[\((n != \ 0)\)\  && \ \((n\  < \ nsoglia)\), \ 
					
        n = nuovagen[p0, p1, n]]; 
	nzeri\  += \ If[n\  == \ 0, \ 1, \ 0]; 
      
		n = 1, 
			\ \ \ {niter}]; 
	N[nzeri/niter]\ \)\)

Costruiamo una funzione che fa tutto ciò in un colpo solo. Anche in questo caso creiamo una funzione compilata.

\(montecarlo = 
      Compile[{p0, p1, {nsoglia, _Integer}, {niter, _Integer}}, 
		
        Module[{n, nzeri = 0, r, p = p0 + p1}, 
			
          Do[
				n = 1; 
				
            While[\((n != \ 0)\) && \((n < nsoglia)\), 
					
              n = Sum[r = Random[]; Which[r < p0, 0, r < p, 1, True, 2], {n}]
                  
					]; 
				If[n == 0, \(nzeri++\), 0], 

            				{niter}]; 
			nzeri/niter\ ], 

        		{{n, _Integer}, {nzeri, _Integer}, {r, _Real}, {p, _Real}, {
            Random[], _Real}}
	]\)

I valori di nsoglia e niter influiscono entrambi sul risultato, ma in modo diverso. Sono soprattuto i casi in cui p0 è piccolo ad essere più critici.

tav = Table[montecarlo[0.1, 0., 10, 100], {10}]

tav = Table[montecarlo[0.1, 0., 10, 1000], {10}]

tav = Table[montecarlo[0.1, 0., 100, 1000], {10}]

tav = Table[montecarlo[0.3, 0., 100, 1000], {10}]

Riferiamoci al caso (Subscript[p, 1]) = 0 e creiamo una tabella con le nostre stime della probabilità di estinzione per valori diversi di (Subscript[p, 0]). In questo esempio poniamo nsoglia=100 e niter=1000, ma questi valori si possono aumentare per ottenere stime più precise.

problis = 
      Table[{x, montecarlo[x, 0., 100, 1000]}, {x, 0., 1., 0.05}]

Possiamo ora rappresentare queste stime in un grafico:

probgraf = ListPlot[problis,  PlotJoined -> True]

Infine possiamo confrontare questo grafico con quello della probabilità esatta che abbiamo ricavato in precedenza

\(Clear[p0, p1, p2]; 
p1 = 0; p2 = 1 - p0 - p1; 

    plim := \ p0/p2\  /; \ p0 < p2; 

    plim\  := \ 1\  /; \ p0 ≥ p2; \)

\(asint = Plot[plim, {p0, 0, 1}, 	PlotStyle -> Hue[1]]\)

Show[probgraf, asint]