Distribuzioni di v.c. continue

Dopo avere introdotto con una simulazione il concetto di funzione di densità nel capitolo La funzione di densità di probabilità, prendiamo ora in esame alcune delle densità di probabilità più note. Iniziamo con la distribuzione continua uniforme.

Preparazione del Notebook

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.stats as stats
import scipy.special as sc
from scipy.integrate import quad
import arviz as az
import preliz as pz

%config InlineBackend.figure_format = 'retina'
RANDOM_SEED = 42
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-darkgrid")
sns.set_theme(palette="colorblind")

Distribuzione uniforme

La distribuzione uniforme è la più sempilce funzione di densità di probabilità. Consideriamo nuovamente l’esperimento con lo spinner che abbiamo introdotto nel capitolo La funzione di densità di probabilità. Simuliamo 20 valori che potrebbero essere ottenuti facendo ruotare lo spinner e li rappresentiamo con un istogramma.

y = rng.uniform(low=0, high=360, size=20)
print(y)

[278.62417748 157.99623831 309.09525117 251.05249046  33.90384524
 351.22404659 274.01029272 282.9831499   46.12090776 162.13893764
 133.48728872 333.63539599 231.79144323 296.19418078 159.62911158
  81.80593984 199.65052333  22.9742122  297.94722192 227.39918368]

plt.figure()
count, bins, ignored = plt.hist(y, bins=36, density=True)
plt.xlabel("Risultato dello spinner")
plt.ylabel("Frequenza relativa");

Sebbene possiamo pensare che sia ugualmente probabile che si verifichi qualsiasi risultato tra 0 e 360, l’istogramma non sembra suggerire questo. Ma lo spinner è stato fatto ruotare solo 20 volte. Proviamo con 100,000 ripetizioni.

plt.figure()
count, bins, ignored = plt.hist(360 * rng.uniform(size=100000), bins=36, density=True)
plt.xlabel("Risultato dello spinner")
plt.ylabel("Frequenza relativa")

Text(0, 0.5, 'Frequenza relativa')

In questo caso, anche se c’è una variazione nelle altezze delle barre (con \(\Delta\) = 10), la forma generale dell’istogramma sembra essere piuttosto piatta, ovvero uniforme, nell’intero intervallo dei valori possibili di \(X\), ovvero \(0 <= X <= 360\). Se potessimo ottenere un numero enorme di risultati dello spinner, il profilo dell’istogramma assumerebbe la forma della funzione di densità uniforme mostratra nella figura seguente.

plt.figure()
x = np.linspace(0, 360, 100)
plt.plot(x, stats.uniform.pdf(x, 0, 360), lw=2, label="uniform pdf")
plt.xlabel("x")
plt.ylabel("p(x)");

La sintassi di preliz è la seguente.

pz.Uniform(lower=0, upper=360).plot_pdf(support=(0, 360));

Quando la variabile casuale \(X\) è continua, come nel caso del risultato della rotazione dello spinner, allora per rappresentare le probabilità usiamo una curva chiamata funzione di densità di probabilità. Poiché la scala dello spinner va da 0 a 360, sappiamo che tutti i risultati possibili devono cadere in questo intervallo, quindi la probabilità che \(X\) assuma un valore nell’intervallo [0, 360] è 1.0. Questa probabilità è rappresentata dall’area totale sotto la funzione di densità della figura precedente tra 0 e 360. Poiché l’area di questo rettangolo è data dall’altezza per la base e la base è uguale a 360, l’altezza di questa curva di densità deve essere 1/360 = 0.00278. L’ordinata della funzione di densità (qui 0.00278 nell’intervallo [0, 360] e 0 altrove) è chiamata densità.

Le probabilità corrispondono alle aree sottese alla curva di densità nell’intervallo di valori \(X\) specificato. Per esempio, nell’esperimento dello spinner possiamo chiederci quale sia la probabilità di ottenere un numero compreso tra 150 e 250, ovvero \(P(150 < X < 250)\). Per trovare la risposta dobbiamo calcolare l’area di un rettangolo. La base è 250 - 150 = 100. L’altezza è 0.00278. Dunque, la probabilità è

100*1/360

0.2777777777777778

Per svolgere questo calcolo i software utilizzano la funzione di ripartizione, \(P(X < x)\). Per trovare l’area in un intervallo è necessario sottrarre due aree. Nel caso presente abbiamo \(P(x < 250) - P(x < 150)\), ovvero:

stats.uniform.cdf(250, 0, 360) - stats.uniform.cdf(150, 0, 360)

0.27777777777777773

La probabilità cercata è rappresentata dal rettangolo indicato nella figura seguente.

plt.figure()
x = np.linspace(0, 360, 1000)
fx = stats.uniform.pdf(x, 0, 360)
plt.plot(x, fx)
plt.fill_between(x, fx, where=(x >= 150) & (x <= 250), color="0.75");

In maniera più formale possiamo dire che la distribuzione continua uniforme è una distribuzione di probabilità continua che assegna lo stesso grado di fiducia a tutti i possibili valori di una variabile definita in un certo intervallo \(S=[a,b]\subset {\mathbb {R}}\). La distribuzione continua uniforme viene indicata con \({\mathcal {U}}(a,b)={\mathcal {U}}([a,b])\). Come intervallo \([a,b]\) viene spesso preso l’intervallo unitario \(I=[0,1]\).

La densità di probabilità di una variabile casuale continua uniforme \({\mathcal {U}}(a,b)\) è

\[ f(x)={\frac {1}{b-a}} \quad \text{su}\; [a, b]. \]

Il suo valore attesto è

\[ \displaystyle E(X)={\frac {1}{2}}(b+a). \]

La sua varianza è

\[ V(X)={\frac {1}{12}}(b-a)^{2}. \]

In Python è possibile manipolare la distribuzione uniforme mediante la funzione uniform del modulo scipy.stats. Di default, la funzione scipy.stats.uniform() è un’istanziazione di \({\mathcal{U}}(0,1)\). Se utilizziamo la funzione pdf() (probability density function) otteniamo l’ordinata della funzione di densità \({\mathcal{U}}(0,1)\) in corrispondenza dei valori \(x\) passati in input. Per esempio, esaminiamo la funzione di densità \({\mathcal{U}}(0,1)\) in corrispondenza di 0.5, 0.8 e 1.2. Per i primi due valori ci aspettiamo di ottenere 1; in corrispondenza di 1.2 ci aspettiamo di ottenere 0, poiché questo valore è al di fuori dell’intervallo \([ 0, 1]\).

stats.uniform.pdf([0.5, 0.8, 1.2])

array([1., 1., 0.])

Con la funzione cdf() (cumulative density function) otteniamo la funzione di ripartizione. Per esempio, per \({\mathcal{U}}(0,1)\) in corrispondenza dei punti 0.5 e 0.8 otteniamo

stats.uniform.cdf([0.5, 0.8])

array([0.5, 0.8])

Usando la funzione di ripartizione è possibile calcolare la probabilità che la variabile casuale continua assuma un valore nell’intervallo specificato. Per esempio, per \({\mathcal{U}}(0,1)\) troviamo \(P(0.5 < x < 0.8)\)

stats.uniform.cdf(0.8) - stats.uniform.cdf(0.5)

0.30000000000000004

I quantili di una funzione di densità (ovvero, il valore della variabile casuale \(X\) in corrispondenza del valore della funzione di ripartizione fornito in input) si ottengono con la funzione ppf() (probability point function). Per esempio, troviamo i quantili di ordine 0.5 e 0.8 di una \({\mathcal {U}}(0,1)\).

stats.uniform.ppf([0.5, 0.8])

array([0.5, 0.8])

Infine, è possibile simulare dei valori casuali della distribuzione \({\mathcal{U}}(0,1)\) usando la funzione rng.uniform(). Se vogliamo 5 valori da una \({\mathcal{U}}(0,1)\), scriviamo:

rng.uniform(size=5)

array([0.37949106, 0.27866256, 0.8162908 , 0.69669817, 0.10985629])

Verifico il valore atteso di 100,000 realizzazioni di \({\mathcal {U}}(0,1)\).

rng.uniform(size=100000).mean()

0.49932833719395314

Verifico la varianza di 100,000 realizzazioni di \({\mathcal {U}}(0,1)\).

rng.uniform(size=100000).var()

0.0832116489056256

1 / 12

0.08333333333333333

Distribuzione esponenziale

Un’altra distribuzione di densità molto semplice è la distribuzione esponenziale. La distribuzione esponenziale viene spesso utilizzata per modellare il tempo trascorso prima che un evento si verifichi (tempo di attesa).

La distribuzione esponenziale è l’unica distribuzione di probabilità continua che possiede la proprietà di assenza di memoria. Ad esempio, ipotizziamo che il tempo necessario affinché un bicchiere da vino si rompa dopo il primo utilizzo segua una distribuzione esponenziale. Supponiamo inoltre che ci sia un bicchiere da vino che non si è rotto dopo 3 anni dal primo utilizzo. L’assenza di memoria significa che la probabilità che questo bicchiere da vino non si rompa nel prossimo anno è la stessa della probabilità che un altro bicchiere da vino nuovo non si rompa nel primo anno di utilizzo.

Chiamiamo \(X\) il tempo di attesa. Sia \(\mu = \mathbb{E}(X)\) il tempo di attesa medio. La funzione di densità esponenziale è

(41)\[ f(x) = \lambda {\rm e}^{-\lambda x}, \quad \text{con} \; \lambda = 1/\mu,\, \lambda > 0,\, x > 0, \]

ovvero

\[ f(x) = \frac{1}{\mu} {\rm e}^{-x/\mu}. \]

La media di una distribuzione esponenziale è

\[ E(X) = \frac{1}{\lambda}. \]

La varianza di una distribuzione esponenziale è

\[ V(X) = \mu = \frac{1}{\lambda^2}. \]

La deviazione standard è dunque uguale alla media:

\[ \sigma_X = \frac{1}{\lambda} = \mu. \]

Ad esempio, il tempo di attesa della pubblicazione del voto di un esame scritto segue una distribuzione esponenziale. Supponiamo che, in questo Corso di Laurea, il tempo di attesa medio per conoscere il risultato di un esame scritto sia di 4 giorni. La funzione esponenziale diventa

\[ f(x) = \frac{1}{4} \exp^{-x/4}. \]

Per disegnare un grafico della funzione esponenziale possiamo usare la funzione scipy.stats.expon(). La densità è data da pdf(x, loc, scale), laddove il parametro loc è 0 e scale è la deviazione standard. Nel caso presente abbiamo

x = np.arange(0, 20, 0.01)
mu = 4
lam = 1 / mu
stdev = 1 / lam
pdf = stats.expon.pdf(x, loc=0, scale=stdev)

plt.figure()
plt.plot(x, pdf)
plt.xlabel("x")
plt.ylabel("f(x)");

Chiediamoci, ad esempio quale sia la probabilità di dovere aspettare non più di un giorno e mezzo per conoscere il voto. La risposta a questa domanda è data dalla funzione di ripartizione in corrispondenza di 1.5, ovvero \(F(1.5) = P(X \leq 1.5)\).

fx = stats.expon.pdf(x, loc=0, scale=stdev)

plt.figure()
plt.plot(x, fx)
plt.fill_between(x, fx, where=(x >= 0) & (x <= 1.5), color="0.75");

Possiamo trovare la risposta usando la funzione cdf():

stats.expon.cdf(1.5, loc=0, scale=stdev) 

0.3127107212090278

Chiediamoci, ad esempio quale sia la probabilità di conoscere il voto dopo avere aspettato un tempo compreso tra 1 e 6 giorni. Dobbiamo trovare l’area sottesa alla funzione di densità nell’intervallo [1, 6]. Usando la fuzione di ripartizione, calcoliamo \(F(6) - F(1) = P(X <= 6) - P(X <= 1)\).

fx = stats.expon.pdf(x, loc=0, scale=stdev)

plt.figure()
plt.plot(x, fx)
plt.fill_between(x, fx, where=(x >= 1) & (x <= 6), color="0.75");

stats.expon.cdf(6, loc=0, scale=stdev) - stats.expon.cdf(1, loc=0, scale=stdev)

0.5556706229229751

Troviamo la probabilità di dovere aspettare almeno 5 giorni e mezzo.

fx = stats.expon.pdf(x, loc=0, scale=stdev)

plt.figure()
plt.plot(x, fx)
plt.fill_between(x, fx, where=(x >= 5.5) & (x <= 21), color="0.75");

La probabilità cercata è data dalla probabilità dell’evento complementare di quello fornito dalla funzione di ripartizione.

1 - stats.expon.cdf(5.5, loc=0, scale=stdev) 

0.25283959580474646

Se la media del tempo di attesa nel Corso di Laurea fosse di 4 giorni, allora circa una volta su 4 lo studente dovrà aspettare almeno 5.5 giorni per conoscere il voto dello scritto.

La figura seguente mostra un istogramma di 1000000 valori casuali estratti dalla distribuzione esponenziale di parametro \(\lambda = 1/4\). All’istogramma è sovrapposta la funzione di densità.

sexpon = stats.expon(loc=0, scale=stdev)
n = 1000000
samps = sexpon.rvs(size=n)

plt.figure()
count, bins, ignored = plt.hist(samps, bins=100, density=True)
plt.plot(x, fx)
plt.xlim([0, 20])
plt.ylabel("Frequenza relativa")
plt.xlabel("Tempo di attesa");

pz.Exponential(lam = 1/4).plot_pdf(support=(0, 20));

Distribuzione Gaussiana

La più importante distribuzione di densità è la Gaussiana. Non c’è un’unica distribuzione gaussiana (o Normale): la distribuzione gaussiana è una famiglia di distribuzioni. Tali distribuzioni sono dette “gaussiane” in onore di Carl Friedrich Gauss (uno dei più grandi matematici della storia il quale, tra le altre cose, scoprì l’utilità di tale funzione di densità per descrivere gli errori di misurazione). Adolphe Quetelet, il padre delle scienze sociali quantitative, fu il primo ad applicare tale funzione di densità alle misurazioni dell’uomo. Karl Pearson usò per primo il termine “distribuzione normale” anche se ammise che questa espressione “ha lo svantaggio di indurre le persone a credere che le altre distribuzioni, in un senso o nell’altro, non siano normali.”

Limite delle distribuzioni binomiali

Iniziamo con un un breve excursus storico. Nel 1733, Abraham de Moivre notò che, aumentando il numero di prove di una distribuzione binomiale, la distribuzione risultante diventava quasi simmetrica e a forma campanulare. Con 10 prove e una probabilità di successo di 0.9, la distribuzione è chiaramente asimmetrica.

n = 10
p = 0.9
r_values = list(range(n + 1))
dist = [stats.binom.pmf(r, n, p) for r in r_values]

plt.figure()
plt.bar(r_values, dist);

Ma se aumentiamo il numero di prove di un fattore di 100 a N = 1000, senza modificare la probabilità di successo di 0.9, la distribuzione assume una forma campanulare quasi simmetrica. Dunque, de Moivre scoprì che, quando N è grande, la funzione gaussiana (che introdurremo qui sotto), nonostante sia la densità di v.c. continue, fornisce una buona approssimazione alla funzione di massa di probabilità binomiale.

n = 1000
p = 0.9
r_values = list(range(n + 1))
dist = [stats.binom.pmf(r, n, p) for r in r_values]

plt.figure()
plt.bar(r_values, dist)
plt.xlim(850, 950);

La distribuzione Normale fu scoperta da Gauss nel 1809. Il Paragrafo successivo illustra come si possa giungere alla Normale mediante una simulazione.

La Normale prodotta con una simulazione

Il libro “Rethinking Statistics” di McElreath [McE20] spiega come sia possibile ottenere la distribuzione normale attraverso una simulazione. Immaginiamo di avere duemila persone che si trovano allineate su una linea di partenza. Quando viene dato il segnale di partenza, ogni persona lancia una moneta e compie un passo avanti o indietro a seconda del risultato del lancio. La lunghezza di ogni passo può variare da 0 a 1 metro. Ogni persona lancia la moneta 16 volte e quindi compie 16 passi.

def randomwalk(n):
    steps = []
    for i in range(n):
        rand = rng.integers(1, 3)
        if rand == 1:
            steps.append(-rng.uniform(0, 1))
        else:
            steps.append(rng.uniform(0, 1))
    walk = np.cumsum(steps)
    return walk

I valori ottenuti da una singola passeggiata casuale consistono nella distanza dall’origine (punto di partenza, indicato come 0) dopo un determinato numero di passi. Questi valori sono rappresentati da numeri. Per esempio:

particularWalk = randomwalk(16)
print(*particularWalk)

0.37531838354697544 1.0584896463887754 1.4534636834009325 1.2937899514171627 1.5447231066244227 1.3373023369291819 1.3254604685050029 0.7344137731826899 1.136894490699304 1.9684525255384093 1.8250625862950245 1.016357414464169 1.231202594793971 2.1991617340098744 2.3950630949221305 2.4013927428252604

Alla conclusione di queste passeggiate casuali (random walk) non possiamo sapere con esattezza dove si troverà ciascuna persona, ma possiamo conoscere con certezza le caratteristiche della distribuzione delle 2000 distanze dall’origine.

for k in range(2000):
    particularWalk = randomwalk(16)
    plt.plot(np.arange(16), particularWalk, linewidth=0.4, alpha=0.1)

plt.plot(np.arange(16), particularWalk, linewidth=0.9)
plt.xlabel("Numero di passi")
plt.ylabel("Distanza dall'origine");

Per esempio, possiamo predire in maniera accurata la proporzione di persone che si sono spostate in avanti oppure all’indietro. Oppure, possiamo predire accuratamente la proporzione di persone che si troveranno ad una certa distanza dalla linea di partenza (es., a 1.5 m dall’origine). Queste predizioni sono possibili perché tali distanze si distribuiscono secondo la legge Normale (la curva rossa nel grafico seguente).

distances = []
for k in range(2000):
    distances.append(randomwalk(16)[15])

avg_distances = np.mean(distances)
std_distances = np.std(distances)

plt.figure()
plt.hist(
    distances, density=True, bins=20, label="distanze", alpha=0.7
)

x_pdf = np.linspace(min(distances), max(distances), 100)
y_pdf = stats.norm.pdf(x_pdf, avg_distances, std_distances)
plt.plot(x_pdf, y_pdf, "C0", lw=2, label="pdf")

plt.xlabel("Distanza percorsa (m)")
plt.ylabel("Densità")
plt.legend();

Questa simulazione dimostra che qualsiasi processo che prevede la somma di un certo numero di valori casuali, tutti provenienti dalla stessa distribuzione, alla fine converge ad una distribuzione normale. Non importa quale sia la forma della distribuzione di partenza, che può essere uniforme come nell’esempio mostrato o di qualsiasi altro tipo. La forma della distribuzione da cui viene realizzato il campionamento determina la velocità di convergenza alla distribuzione normale. In alcuni casi, la convergenza può essere lenta, mentre in altri casi, come nell’esempio mostrato, la convergenza può essere molto rapida. È possibile osservare lo stesso principio tramite la Galton box.

Dal punto di vista formale, possiamo definire una variabile casuale continua \(Y\) come avente una distribuzione normale se la sua densità di probabilità è distribuita secondo la seguente equazione

(42)\[ f(y; \mu, \sigma) = {1 \over {\sigma\sqrt{2\pi} }} \exp \left\{-\frac{(y - \mu)^2}{2 \sigma^2} \right\}, \]

dove \(\mu \in \mathbb{R}\) e \(\sigma > 0\) sono i parametri della distribuzione.

La densità normale è unimodale e simmetrica con una caratteristica forma a campana e con il punto di massima densità in corrispondenza di \(\mu\).

Il significato dei parametri \(\mu\) e \(\sigma\) che appaiono nell’eq. (42) viene chiarito dalla dimostrazione che

\[ \mathbb{E}(Y) = \mu, \qquad \mathbb{V}(Y) = \sigma^2. \]

La rappresentazione grafica di quattro densità Normali con medie -1, -0.5, 0, 1 e con deviazioni standard 0.25, 0.5, 1 e 2 è fornita nella figura seguente.

x = np.arange(-5, 6, 0.001)

mus = [-1.0, -0.5, 0.0, 1.0]
sigmas = [0.25, 0.5, 1, 2]

plt.figure()

for mu, sigma in zip(mus, sigmas):
    pdf = stats.norm.pdf(x, mu, sigma)
    plt.plot(x, pdf, label=r"$\mu$ = {}, $\sigma$ = {}".format(mu, sigma))
plt.xlabel("x")
plt.ylabel("f(x)")
plt.legend(loc=1);

Concentrazione

È istruttivo osservare il grado di concentrazione della distribuzione Normale attorno alla media:

\[\begin{split} \begin{align} P(\mu - \sigma < Y < \mu + \sigma) &= P (-1 < Z < 1) \simeq 0.683, \notag\\ P(\mu - 2\sigma < Y < \mu + 2\sigma) &= P (-2 < Z < 2) \simeq 0.956, \notag\\ P(\mu - 3\sigma < Y < \mu + 3\sigma) &= P (-3 < Z < 3) \simeq 0.997. \notag \end{align} \end{split}\]

Si noti come un dato la cui distanza dalla media è superiore a 3 volte la deviazione standard presenti un carattere di eccezionalità perché meno del 0.3% dei dati della distribuzione Normale presentano questa caratteristica.

Per indicare la distribuzione Normale si usa la notazione \(\mathcal{N}(\mu, \sigma)\).

Funzione di ripartizione

Il valore della funzione di ripartizione di \(Y\) nel punto \(y\) è l’area sottesa alla curva di densità \(f(y)\) nella semiretta \((-\infty, y]\). Non esiste alcuna funzione elementare per la funzione di ripartizione

(43)\[ F(y) = \int_{-\infty}^y {1 \over {\sigma\sqrt{2\pi} }} \exp \left\{-\frac{(y - \mu)^2}{2\sigma^2} \right\} dy, \]

pertanto le probabilità \(P(Y < y)\) vengono calcolate mediante integrazione numerica approssimata. I valori della funzione di ripartizione di una variabile casuale Normale sono dunque forniti da un software.

Esaminiamo le funzioni fornite da scipy.stats per la densità Normale. Il metodo .norm.rvs() produce un valore casuale estratto dalla distribuzione Normale specificata. Per esempio, un valore casuale dalla \(\mathcal{N}(\mu = 100, \sigma = 15)\) è:

stats.norm.rvs(loc=100, scale=15)

102.31135123904185

Estraiamo ora 10 valori a caso dalla \(\mathcal{N}(100, 15)\):

qi = stats.norm.rvs(loc=100, scale=15, size=10)
print(*qi)

66.64822051568693 88.90181634038653 99.51168016916533 111.35374005701324 106.53766909647265 112.84448406807662 103.65149109665887 89.79044973299285 112.05607385093361 101.43519453618879

Per trovare la probabilità che un’osservazione estratta a caso dalla \(\mathcal{N}(100, 15)\) abbia un valore minore o uguale a, diciamo, 115, troviamo il valore della funzione di ripartizione (o funzione cumulativa di densità) nel punto 115.

stats.norm.cdf(115, 100, 15)

0.8413447460685429

Questa è l’area sottesa alla funzione di densità nell’intervallo \([-\infty, 115]\), come indicato nella figura seguente.

mu = 100
sigma = 15
x = np.linspace(mu - 3 * sigma, mu + 3 * sigma, 10000)
fx = stats.norm.pdf(x, mu, sigma)

plt.figure()
plt.plot(x, fx)
_ = plt.fill_between(x, fx, where=x <= 115, color="0.75")

Solo per fare un esempio, qui di seguito fornisco il codice Python per calcolare l’integrale che stiamo discutendo per mezzo della funzione quad della libreria SciPy:

def gaussian(x, mu, sig):
    return (
        1.0 / (np.sqrt(2.0 * np.pi) * sig) * np.exp(-np.power((x - mu) / sig, 2.0) / 2)
    )

mu = 100
sigma = 15
result, error = quad(gaussian, -1000, 115, args=(mu, sigma))
print("Il risultato è", result, "con errore", error)

Il risultato è 0.841344746068543 con errore 4.0191205799653545e-10

Il risultato replica quello prodotto da .norm.cdf().

Per trovare la proporzione di persone nella popolazione che hanno un QI maggiore di 2 deviazioni standard dalla media consideriamo l’evento complementare:

1 - stats.norm.cdf(130, 100, 15)

0.02275013194817921

plt.figure()
plt.plot(x, fx)
plt.fill_between(x, fx, where=x >= 130, color="0.75");

In maniera equivalente, possiamo usare la Survival Function:

stats.norm.sf(130, 100, 15)

0.022750131948179195

La funzione ppf restituisce il quantile della Normale. Ad esempio:

stats.norm.ppf(1 - 0.022750131948179195, 100, 15)

130.0

Distribuzione Normale standard

La distribuzione Normale di parametri \(\mu = 0\) e \(\sigma = 1\) viene detta distribuzione Normale standard. La famiglia Normale è l’insieme avente come elementi tutte le distribuzioni Normali con parametri \(\mu\) e \(\sigma\) diversi. Tutte le distribuzioni Normali si ottengono dalla Normale standard mediante una trasformazione lineare: se \(Y \sim \mathcal{N}(\mu_Y, \sigma_Y)\) allora

\[ X = a + b Y \sim \mathcal{N}(\mu_X = a+b \mu_Y, \sigma_X = \left|b\right|\sigma_Y). \]

L’area sottesa alla curva di densità di \(\mathcal{N}(\mu, \sigma)\) nella semiretta \((-\infty, y]\) è uguale all’area sottesa alla densità Normale standard nella semiretta \((-\infty, z]\), in cui \(z = (y -\mu_Y )/\sigma_Y\) è il punteggio standard di \(Y\). Per la simmetria della distribuzione, l’area sottesa nella semiretta \([1, \infty)\) è uguale all’area sottesa nella semiretta \((-\infty, 1]\) e quest’ultima coincide con \(F(-1)\). Analogamente, l’area sottesa nell’intervallo \([y_a, y_b]\), con \(y_a < y_b\), è pari a \(F(z_b) - F(z_a)\), dove \(z_a\) e \(z_b\) sono i punteggi standard di \(y_a\) e \(y_b\).

Si ha anche il problema inverso rispetto a quello del calcolo delle aree: dato un numero \(0 \leq p \leq 1\), il problema è quello di determinare un numero \(z \in \mathbb{R}\) tale che \(P(Z < z) = p\). Il valore \(z\) cercato è detto quantile di ordine \(p\) della Normale standard e può essere trovato mediante un software.

Supponiamo che l’altezza degli individui adulti segua la distribuzione Normale di media \(\mu = 1.7\) m e deviazione standard \(\sigma = 0.1\) m. Vogliamo sapere la proporzione di individui adulti con un’altezza compresa tra \(1.7\) e \(1.8\) m.

Il problema ci chiede di trovare l’area sottesa alla distribuzione \(\mathcal{N}(\mu = 1.7, \sigma = 0.1)\) nell’intervallo \([1.7, 1.8]\):

stats.norm.cdf(1.8, 1.7, 0.1) - stats.norm.cdf(1.7, 1.7, 0.1)

0.34134474606854315

mu = 1.7
sigma = 0.1
x = np.linspace(mu - 3 * sigma, mu + 3 * sigma, 10000)
fx = stats.norm.pdf(x, mu, sigma)
plt.figure()
plt.plot(x, fx)
plt.fill_between(x, fx, where=(x >= 1.7) & (x <= 1.8), color="0.75");

In maniera equivalente, possiamo standardizzare i valori che delimitano l’intervallo considerato e utilizzare la funzione di ripartizione della normale standardizzata. I limiti inferiore e superiore dell’intervallo sono

\[ z_{\text{inf}} = \frac{1.7 - 1.7}{0.1} = 0, \quad z_{\text{sup}} = \frac{1.8 - 1.7}{0.1} = 1.0, \]

quindi otteniamo

stats.norm.cdf(1.0, 0, 1) - stats.norm.cdf(0, 0, 1)

0.3413447460685429

Il modo più semplice per risolvere questo problema resta comunque quello di rendersi conto che la probabilità richiesta non è altro che la metà dell’area sottesa dalle distribuzioni Normali nell’intervallo \([\mu - \sigma, \mu + \sigma]\), ovvero \(0.683/2\).

Si noti che, usando preliz, è molto semplice ottenere un grafico della funzione gaussiana. Per esempio, nel caso della normale standard abbiamo:

pz.Normal(0, 1).plot_pdf(pointinterval=True);

Distribuzione Chi-quadrato

Dalla Normale deriva la distribuzione \(\chi^2\). La distribuzione \(\chi^2_{~k}\) con \(k\) gradi di libertà descrive la variabile casuale

\[ Z_1^2 + Z_2^2 + \dots + Z_k^2, \]

dove \(Z_1, Z_2, \dots, Z_k\) sono variabili casuali i.i.d. che seguono la distribuzione Normale standard \(\mathcal{N}(0, 1)\). La variabile casuale chi-quadrato dipende dal parametro intero positivo \(\nu = k\) che ne identifica il numero di gradi di libertà. La densità di probabilità di \(\chi^2_{~\nu}\) è

\[ f(x) = C_{\nu} x^{\nu/2-1} \exp (-x/2), \qquad \text{se } x > 0, \]

dove \(C_{\nu}\) è una costante positiva.

La figura seguente mostra alcune distribuzioni Chi-quadrato variando il parametro \(\nu\).

x = np.arange(0, 40, 0.1)

nus = [2, 4, 8, 16]
plt.figure()
for nu in nus:
    pdf = stats.chi2.pdf(x, nu)
    plt.plot(x, pdf, label=r"$\nu$ = {}".format(nu))
plt.xlabel("x")
plt.ylabel("f(x)")
plt.legend(loc=1);

Proprietà

• La distribuzione di densità \(\chi^2_{~\nu}\) è asimmetrica.

• Il valore atteso di una variabile \(\chi^2_{~\nu}\) è uguale a \(\nu\).

• La varianza di una variabile \(\chi^2_{~\nu}\) è uguale a \(2\nu\).

• Per \(k \rightarrow \infty\), la \(\chi^2_{~\nu} \rightarrow \mathcal{N}\).

• Se \(X\) e \(Y\) sono due variabili casuali chi-quadrato indipendenti con \(\nu_1\) e \(\nu_2\) gradi di libertà, ne segue che \(X + Y \sim \chi^2_m\), con \(m = \nu_1 + \nu_2\). Tale principio si estende a qualunque numero finito di variabili casuali chi-quadrato indipendenti.

Per fare un esempio, consideriamo la v.c. \(\chi^2_5\).

# Set the degrees of freedom
df = 5

# Create a chi-square distribution object
chi2_dist = stats.chi2(df)

# Generate x values for the plot
x = np.linspace(0, 20, 200)

# Calculate the probability density function (PDF) of the chi-square distribution for x values
pdf = chi2_dist.pdf(x)

# Plot the PDF
plt.figure()
plt.plot(x, pdf)
plt.title('Chi-Square Distribution (df=5)')
plt.xlabel('x')
plt.ylabel('PDF');

Generiamo 1000000 valori da questa distribuzione.

x = chi2_dist.rvs(1000000)
x[0:20]

array([ 1.45264069,  2.18922686,  5.85814879,  7.56103994,  5.15554681,
        7.50144358,  5.18620004,  3.41722509,  2.55049895,  3.06276035,
        3.14751315,  2.74953353, 10.94120946,  0.35083943,  4.91169935,
        6.32132976,  4.08343197,  1.44730045,  4.24364635,  9.07363732])

Calcoliamo la media di questi valori.

np.mean(x)

4.998638907259898

Calcolo la varianza.

np.var(x, ddof=0)

9.988859276091533

Usando preliz, possiamo produrre un grafico della funzione \(\chi^2_5\) nel modo seguente:

pz.ChiSquared(5).plot_pdf(pointinterval=True);

Distribuzione \(t\) di Student

Dalle distribuzioni Normale e Chi-quadrato deriva un’altra distribuzione molto nota, la \(t\) di Student. Se \(Z \sim \mathcal{N}\) e \(W \sim \chi^2_{~\nu}\) sono due variabili casuali indipendenti, allora il rapporto

(44)\[ T = \frac{Z}{\Big( \frac{W}{\nu}\Big)^{\frac{1}{2}}} \]

definisce la distribuzione \(t\) di Student con \(\nu\) gradi di libertà. Si usa scrivere \(T \sim t_{\nu}\). L’andamento della distribuzione \(t\) di Student è simile a quello della distribuzione Normale, ma ha una dispersione maggiore (ha le code più pesanti di una Normale, ovvero ha una varianza maggiore di 1).

La seguente mostra alcune distribuzioni \(t\) di Student variando il parametro \(\nu\).

x = np.arange(-5, 5, 0.1)

nus = [1, 2, 5, 30]

plt.figure()
for nu in nus:
    pdf = stats.t.pdf(x, nu)
    plt.plot(x, pdf, label=r"$\nu$ = {}".format(nu))
plt.plot(x, stats.norm.pdf(x, 0, 1), label="N(μ = 0, σ = 1)")
plt.xlabel("x")
plt.ylabel("f(x)")
plt.legend(loc=1);

Proprietà

La variabile casuale \(t\) di Student soddisfa le seguenti proprietà:

1. Per \(\nu \rightarrow \infty\), \(t_{\nu}\) tende alla normale standard \(\mathcal{N}(0, 1)\).

2. La densità della \(t_{\nu}\) è una funzione simmetrica con valore atteso nullo.

3. Per \(\nu > 2\), la varianza della \(t_{\nu}\) vale \(\nu/(\nu - 2)\); pertanto è sempre maggiore di 1 e tende a 1 per \(\nu \rightarrow \infty\).

Per esempio, calcoliamo il valore della funzione di ripartizione di ordine 0.025 nel caso di una \(t_{30}\).

stats.t.ppf(0.025, 30)

-2.042272456301238

Aumentiamo i gradi di libertà: \(\nu\) = 1000.

stats.t.ppf(0.025, 1000)

-1.9623390808264078

Questo valore è quasi identico a quello della Normale stanardizzata.

stats.norm.ppf(0.025, 0, 1)

-1.9599639845400545

La ragione per cui il quantile della distribuzione \(t\) con \(\nu=30\) è maggiore (in valore assoluto) del quantile omotetico della distribuzione Normale Standard è che la distribuzione \(t\) ha una varianza maggiore rispetto alla distribuzione Normale Standard.

Usando preliz, la T di Student si visualizza nel modo seguente.

pz.StudentT(nu = 5, mu = 0, sigma = 1).plot_pdf(pointinterval=True);

Funzione beta di Eulero

La funzione beta di Eulero è una funzione matematica, non una densità di probabilità. La menzioniamo qui perché viene utilizzata nella distribuzione Beta. La funzione beta di Eulero si può scrivere in molti modi diversi; per i nostri scopi la presentiamo così:

(45)\[ B(\alpha, \beta) = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha + \beta)}\,, \]

dove \(\Gamma(x)\) è la funzione Gamma, ovvero il fattoriale discendente, cioè

\[ (x-1)(x-2)\ldots (x-n+1)\notag\,. \]

Per esempio, posti \(\alpha = 3\) e \(\beta = 9\), la funzione beta assume il valore

alpha = 3
beta = 9
sc.beta(alpha, beta)

0.00202020202020202

Lo stesso risultato si ottiene con

((2) * (8 * 7 * 6 * 5 * 4 * 3 * 2)) / (11 * 10 * 9 * 8 * 7 * 6 * 5 * 4 * 3 * 2)

0.00202020202020202

ovvero

alpha = 3
beta = 9
sc.gamma(alpha) * sc.gamma(beta) / sc.gamma(alpha + beta)

0.00202020202020202

Distribuzione Beta

La distribuzione Beta rappresenta un importante strumento nell’ambito delle probabilità, essendo una distribuzione di densità specificamente concepita per descrivere le variabilità delle probabilità. Questa distribuzione è particolarmente utile per modellare fenomeni espressi in percentuali o proporzioni, con l’importante caratteristica di essere definita esclusivamente nell’intervallo \((0, 1)\). In altre parole, essa contempla valori compresi strettamente tra 0 e 1, escludendo questi estremi.

Definizione Formale

Consideriamo una variabile casuale \(\pi\), la quale può assumere qualunque valore nell’intervallo aperto \((0, 1)\). Se diciamo che \(\pi\) segue una distribuzione Beta con parametri \(\alpha\) e \(\beta\) (indicato come \(\pi \sim \text{Beta}(\alpha, \beta)\)), intendiamo che la sua funzione di densità è descritta dalla seguente formula:

\[ \text{Beta}(\pi \mid \alpha, \beta) = \frac{1}{B(\alpha, \beta)}\pi^{\alpha-1} (1-\pi)^{\beta-1} = \frac{\Gamma(\alpha+ \beta)}{\Gamma(\alpha)\Gamma(\beta)}\pi^{\alpha-1} (1-\pi)^{\beta-1} \quad \text{per } \pi \in (0, 1)\,, \]

dove \(B(\alpha, \beta)\) è la funzione beta di Eulero.

I Parametri \(\alpha\) e \(\beta\)

I parametri \(\alpha\) e \(\beta\) giocano un ruolo cruciale nella distribuzione Beta, influenzando direttamente la sua forma e il suo comportamento. È essenziale che entrambi questi parametri siano positivi.

Intuizione e Collegamento con la Distribuzione Binomiale

La distribuzione Beta può essere meglio compresa quando la si osserva in relazione con la distribuzione binomiale. Mentre la distribuzione binomiale modella il numero di successi in una serie di prove, la distribuzione Beta si focalizza sulla probabilità di successo in queste prove.

Nel contesto della distribuzione binomiale, la probabilità di successo è un parametro fisso; nella distribuzione Beta, questa probabilità diventa una variabile aleatoria.

Interpretazione dei Parametri \(\alpha\) e \(\beta\)

I parametri \(\alpha\) e \(\beta\) possono essere interpretati come rappresentanti il numero di successi e insuccessi, rispettivamente. Questa interpretazione è analoga ai termini \(n\) e \(n-x\) nella distribuzione binomiale.

La scelta di \(\alpha\) e \(\beta\) dipende dall’aspettativa iniziale della probabilità di successo:

• Se si presume un’alta probabilità di successo (ad esempio, 90%), si potrebbe scegliere \(\alpha = 90\) e \(\beta = 10\).

• Al contrario, per una bassa aspettativa di successo, si potrebbe impostare \(\alpha = 10\) e \(\beta = 90\).

Un aumento di \(\alpha\) (successi) sposta la distribuzione verso destra, mentre un aumento di \(\beta\) (insuccessi) la sposta verso sinistra. Inoltre, se sia \(\alpha\) sia \(\beta\) aumentano, la distribuzione diventa più stretta, indicando una maggiore certezza.

Questa interpretazione consente di utilizzare la distribuzione Beta per esprimere le nostre credenze a priori riguardo a una sequenza di prove di Bernoulli, dove il rapporto tra successi e tentativi totali è dato da:

\[ \frac{\text{Numero di successi}}{\text{Numero di successi} + \text{Numero di insuccessi}} = \frac{\alpha}{\alpha + \beta}\notag\,. \]

Al variare di \(\alpha\) e \(\beta\) si ottengono molte distribuzioni di forma diversa; un’illustrazione è fornita dalla seguente GIF animata.

La figura seguente mostra la distribuzione \(Beta(x \mid \alpha, \beta)\) per \(\alpha\) = 0.5, 5.0, 1.0, 2.0, 2.0 e \(\beta\) = 5, 1.0, 3.0, 2.0, 5.0.

x = np.linspace(0, 1, 200)
alphas = [0.5, 5.0, 1.0, 2.0, 2.0]
betas = [0.5, 1.0, 3.0, 2.0, 5.0]

plt.figure()
for a, b in zip(alphas, betas):
    pdf = stats.beta.pdf(x, a, b)
    plt.plot(x, pdf, label=r"$\alpha$ = {}, $\beta$ = {}".format(a, b))
plt.xlabel("x", fontsize=12)
plt.ylabel("f(x)", fontsize=12)
plt.ylim(0, 4.5)
plt.legend(loc=9);

Costante di normalizzazione

Il rapporto \(\frac{1}{B(\alpha, \beta)} = \frac{\Gamma(\alpha+b)}{\Gamma(\alpha)\Gamma(\beta)}\) è una costante di normalizzazione:

\[ \int_0^1 \pi^{\alpha-1} (1-\pi)^{\beta-1} = \frac{\Gamma(\alpha+b)}{\Gamma(\alpha)\Gamma(\beta)}\,. \]

Ad esempio, con \(\alpha = 3\) e \(\beta = 9\) abbiamo

def integrand(p, a, b):
    return p ** (a - 1) * (1 - p) ** (b - 1)


a = 3
b = 9
result, error = quad(integrand, 0, 1, args=(a, b))
print(result)

0.00202020202020202

ovvero

import math

a = 3
b = 9

result = math.gamma(a) * math.gamma(b) / math.gamma(a + b)
print(result)

0.00202020202020202

ovvero

sc.beta(a, b)

0.00202020202020202

Proprietà

Il valore atteso, la moda e la varianza di una distribuzione Beta sono dati dalle seguenti equazioni:

(46)\[ \mathbb{E}(\pi) = \frac{\alpha}{\alpha+\beta}\,, \]

(47)\[ Mo(\pi) = \frac{\alpha-1}{\alpha+\beta-2}\,, \]

(48)\[ \mathbb{V}(\pi) = \frac{\alpha \beta}{(\alpha+\beta)^2 (\alpha+\beta+1)}\,. \]

La funzione beta_mean_mode_variance() ci restituisce la media, moda e varianza della distribuzione Beta.

def beta_mean_mode_variance(alpha, beta):
    mean = alpha / (alpha + beta)
    mode = (alpha - 1) / (alpha + beta - 2)
    variance = alpha * beta / ((alpha + beta) ** 2 * (alpha + beta + 1))
    return mean, mode, variance

Per esempio

alpha = 3
beta = 9
mean, mode, variance = beta_mean_mode_variance(alpha, beta)
print(f"Mean: {mean}, Mode: {mode}, Variance: {variance}")

Mean: 0.25, Mode: 0.2, Variance: 0.014423076923076924

Usando preliz otteniamo:

pz.Beta(alpha = 2, beta = 5).plot_pdf(pointinterval=True);

Distribuzione a priori coniugata

La distribuzione Beta è la prior coniugata per le distribuzioni Bernoulli, binomiale, binomiale negativa e geometrica (che sono le distribuzioni che coinvolgono successi e fallimenti) nell’inferenza bayesiana.

Calcolare una posterior utilizzando una prior coniugata è molto conveniente, poiché si possono evitare costose computazioni numeriche coinvolte nell’inferenza bayesiana.

Ad esempio, la distribuzione Beta è una prior coniugata per la distribuzione binomiale. Se scegliamo di utilizzare la distribuzione Beta Beta(α, β) come prior, durante la fase di modellazione, sappiamo già che la posterior sarà anch’essa una distribuzione Beta. Pertanto, dopo avere osservato i dati, è possibile calcolare la posterior semplicemente sommando il numero di successi (x) e insuccessi (n-x) ai parametri esistenti α e β rispettivamente, invece di moltiplicare la verosimiglianza per la distribuzione prior. La posterior diventa una distribuzione Beta con parametri (x+α, n-x+β).

Warning

Attenzione alle parole: in questo contesto, il termine “beta” viene utilizzato con tre significati diversi:

• la distribuzione di densità Beta,

• la funzione matematica beta,

• il parametro \(\beta\).

Distribuzione di Cauchy

La distribuzione di Cauchy è un caso speciale della distribuzione di \(t\) di Student con 1 grado di libertà. È definita da una densità di probabilità che corrisponde alla seguente funzione, dipendente da due parametri \(\alpha\) e \(\beta\),

(49)\[ f(x \mid \alpha, \beta) = \frac{1}{\pi \beta \left[1 + \left( \frac{x - \alpha}{\beta} \right)^2\right]}. \]

Il grafico mostra alcune distribuzioni di Cauchy con \(\alpha\) = 0., 0., 0., -2.0 e \(\beta\) = .5, 1., 2., 1.0.

x = np.linspace(-5, 5, 500)
alphas = [0.0, 0.0, 0.0, -2.0]
betas = [0.5, 1.0, 2.0, 1.0]

plt.figure()
for a, b in zip(alphas, betas):
    pdf = stats.cauchy.pdf(x, loc=a, scale=b)
    plt.plot(x, pdf, label=r"$\alpha$ = {}, $\beta$ = {}".format(a, b))
plt.xlabel("x", fontsize=12)
plt.ylabel("f(x)", fontsize=12)
plt.legend(loc=1);

Usando preliz otteniamo:

pz.Cauchy(alpha = 0, beta = 1).plot_pdf(support=(-5,5));

Distribuzione Gamma

La distribuzione gamma è una distribuzione di probabilità continua che gioca un ruolo chiave nella modellazione del tempo di attesa per l’occorrenza di un certo numero di eventi indipendenti e rari. Essa è caratterizzata da due parametri principali: \( \alpha \) e \( \beta \), noti come parametro di forma e parametro di scala, rispettivamente.

Parametro \( \alpha \) - parametro di forma

Il parametro di forma, \( \alpha \), determina la forma generale della curva della distribuzione:

• Se \( \alpha = 1 \), la distribuzione gamma si riduce a una distribuzione esponenziale.

• Se \( \alpha > 1 \), la distribuzione presenta un picco attorno a \( (\alpha - 1) \cdot \beta \).

• Se \( \alpha < 1 \), la distribuzione è inclinata verso destra, mostrando una coda lunga che si estende verso valori più bassi.

In termini interpretativi, \( \alpha \) rappresenta il numero di eventi che si stanno aspettando. Ad esempio, potrebbe rappresentare il numero di ricordi vividi che ci si aspetta di esperire in un certo periodo di tempo.

Parametro \( \beta \) - parametro di scala

Il parametro di scala, \( \beta \), controlla la “larghezza” della distribuzione:

• Un valore più grande di \( \beta \) produce una curva più piatta, indicando una maggiore variabilità nel tempo di attesa.

• Un valore più piccolo di \( \beta \) rende la curva più appuntita, indicando una minore variabilità.

Nel contesto del tempo di attesa, \( \beta \) agisce come una scala temporale, con un valore più grande che indica un periodo di tempo più lungo tra gli eventi, e un valore più piccolo che indica un periodo di tempo più breve.

Formula della funzione di densità di probabilità

La formula matematica per la funzione di densità di probabilità (PDF) della distribuzione gamma è:

\[ f(x|\alpha, \theta) = \frac{x^{\alpha-1}e^{-\frac{x}{\theta}}}{\theta^\alpha \Gamma(\alpha)}, \]

dove,

• \(x\) è la variabile casuale continua, con \(x > 0\).

• \(\theta = \frac{1}{\beta}\).

• \(\Gamma(\alpha)\) è la funzione Gamma di Eulero, che estende la nozione di fattoriale ai numeri reali e complessi. Per numeri interi \(n\), si ha che \(\Gamma(n) = (n-1)!\), ma per argomenti generali \(\alpha\), la funzione Gamma è definita come \(\Gamma(\alpha) = \int_0^\infty x^{\alpha-1}e^{-x}dx\).

Le espressioni per media e varianza della distribuzione Gamma sono le seguenti:

• La media (\(\mu\)) della distribuzione Gamma è \(\mu = \alpha / \beta\), o equivalentemente \(\mu = \alpha \theta\), usando il parametro di scala.

• La varianza (\(\sigma^2\)) della distribuzione Gamma è \(\sigma^2 = \alpha / \beta^2\), o equivalentemente \(\sigma^2 = \alpha \theta^2\), adottando il parametro di scala.

Questo chiarisce la relazione tra i parametri \(\alpha\) (forma) e \(\beta\) (tasso) o \(\theta\) (scala), e come influenzano la distribuzione Gamma.

Per esempio, qui è riportata la distribuzione Gamma di parametri \(\alpha\) = 3 e \(\beta\) = 5/3.

alpha = 3
beta = 5/3

mean = alpha / beta
print(mean)

1.7999999999999998

# Standard deviation = sqrt(alpha / beta^2)

sigma = np.sqrt(alpha / beta**2)
print(sigma)

1.0392304845413263

Creiamo dunque una distribuzione Gamma con una media di tempo d’attesa di 1.8 giorni e deviazione standard pari a 1.04.

pz.Gamma(mu=1.80, sigma=1.04).plot_pdf(pointinterval=True);

pz.Gamma(alpha=3, beta = 5/3).plot_pdf(pointinterval=True);

Distribuzione log-normale

Sia \(y\) una variabile casuale avente distribuzione normale con media \(\mu\) e varianza \(\sigma^2\). Definiamo poi una nuova variabile casuale \(x\) attraverso la relazione

\[ x = e^y \quad \Longleftrightarrow \quad y = \log x. \]

Il dominio di definizione della \(x\) è il semiasse \(x > 0\) e la densità di probabilità \(f(x)\) è data da

(50)\[ f(x) = \frac{1}{\sigma \sqrt{2 \pi}} \frac{1}{x} \exp \left\{-\frac{(\log x - \mu)^2}{2 \sigma^2} \right\}. \]

Questa funzione di densità è chiamata log-normale.

Il valore atteso e la varianza di una distribuzione log-normale sono dati dalle seguenti equazioni:

\[ \mathbb{E}(x) = \exp \left\{\mu + \frac{\sigma^2}{2} \right\}. \]

\[ \mathbb{V}(x) = \exp \left\{2 \mu + \sigma^2 \right\} \left(\exp \left\{\sigma^2 \right\} -1\right). \]

Si può dimostrare che il prodotto di variabili casuali log-normali ed indipendenti segue una distribuzione log-normale.

La figura mostra tre distribuzioni log-normali con \(\mu\) = 0.0 e \(\sigma\) = .25, .5, 1.0.

x = np.linspace(0, 3, 100)
mus = [0.0, 0.0, 0.0]
sigmas = [0.25, 0.5, 1.0]
plt.figure()
for mu, sigma in zip(mus, sigmas):
    pdf = stats.lognorm.pdf(x, sigma, scale=np.exp(mu))
    plt.plot(x, pdf, label=r"$\mu$ = {}, $\sigma$ = {}".format(mu, sigma))
plt.xlabel("x", fontsize=12)
plt.ylabel("f(x)", fontsize=12)
plt.legend(loc=1);

Usando preliz otteniamo:

pz.LogNormal(mu=0, sigma = 0.5).plot_pdf(pointinterval=True);

Commenti e considerazioni finali

È possibile estrarre un campione casuale da diverse distribuzioni utilizzando il generatore di numeri casuali di NumPy. Supponendo che rng sia già stato definito con un seme specifico, possiamo procedere nel modo seguente. Per la distribuzione normale abbiamo:

media, deviazione_standard = 0, 1  # Media e deviazione standard
campione_normale = rng.normal(media, deviazione_standard, size = 100)

laddove size specifica, in questo caso, l’estrazione di 100 valori dalla distribuzione gaussiana. Si procede in maniera simile per altre distribuzioni. Ad esempio

Distribuzione Uniforme:

a, b = 0, 10  # Intervallo
campione_uniforme = rng.uniform(a, b)

Distribuzione t di Student:

gradi_libertà = 10  # Gradi di libertà
campione_t = rng.standard_t(gradi_libertà)

Distribuzione Beta:

alpha, beta = 2, 5  # Parametri
campione_beta = rng.beta(alpha, beta)

Distribuzione Gamma:

forma, scala = 2, 1  # Parametri
campione_gamma = rng.gamma(forma, scala)

Per trovare i valori della funzione di densità di probabilità (PDF) è possibile usare la libreria SciPy. Per esempio, per la normale abbiamo

media, deviazione_standard = 0, 1
x = np.linspace(media - 4*deviazione_standard, media + 4*deviazione_standard, 100)
pdf_normale = stats.norm.pdf(x, loc=media, scale=deviazione_standard)

Altri esempi sono i seguenti.

Distribuzione Uniforme:

a, b = 0, 10
x = np.linspace(a, b, 100)
pdf_uniforme = uniform.pdf(x, loc=a, scale=b-a)

Distribuzione t di Student:

gradi_libertà = 10
x = np.linspace(-5, 5, 100)
pdf_t = stats.t.pdf(x, df=gradi_libertà)

Distribuzione Beta:

alpha, beta_param = 2, 5
x = np.linspace(0, 1, 100)
pdf_beta = stats.beta.pdf(x, alpha, beta_param)

Distribuzione Gamma:

forma, scala = 2, 1
x = np.linspace(0, 10, 100)
pdf_gamma = stats.gamma.pdf(x, a=forma, scale=scala)

La funzione di ripartizione si calcola nel modo seguente per la normale.

media, deviazione_standard = 0, 1
probabilità = 0.5
quantile_normale = stats.norm.ppf(probabilità, loc=media, scale=deviazione_standard)

Altri esempio sono i seguenti.

Distribuzione Uniforme:

a, b = 0, 10
probabilità = 0.5
quantile_uniforme = stats.uniform.ppf(probabilità, loc=a, scale=b-a)

Distribuzione t di Student:

gradi_libertà = 10
probabilità = 0.5
quantile_t = stats.t.ppf(probabilità, df=gradi_libertà)

Distribuzione Beta:

alpha, beta_param = 2, 5
probabilità = 0.5
quantile_beta = stats.beta.ppf(probabilità, alpha, beta_param)

Distribuzione Gamma:

forma, scala = 2, 1
probabilità = 0.5
quantile_gamma = stats.gamma.ppf(probabilità, a=forma, scale=scala)

Infine, il problema opposto, ovvero dato il quantile trovare la probabilità nella coda sinistra della distribuzione, si risolve nel modo seguente. Consideriamo la distribuzione normale.

media, deviazione_standard = 0, 1
quantile = 0
probabilità_normale = stats.norm.cdf(quantile, loc=media, scale=deviazione_standard)

Altri esempi sono i seguenti.

Distribuzione Uniforme:

a, b = 0, 10
quantile = 5
probabilità_uniforme = stats.uniform.cdf(quantile, loc=a, scale=b-a)

Distribuzione t di Student:

gradi_libertà = 10
quantile = 0
probabilità_t = stats.t.cdf(quantile, df=gradi_libertà)

Distribuzione Beta:

alpha, beta_param = 2, 5
quantile = 0.5
probabilità_beta = stats.beta.cdf(quantile, alpha, beta_param)

Distribuzione Gamma:

forma, scala = 2, 1
quantile = 2
probabilità_gamma = stats.gamma.cdf(quantile, a=forma, scale=scala)

%load_ext watermark
%watermark -n -u -v -iv -w

Last updated: Fri Feb 02 2024

Python implementation: CPython
Python version       : 3.11.7
IPython version      : 8.19.0

scipy     : 1.11.4
pandas    : 2.1.4
numpy     : 1.26.2
seaborn   : 0.13.0
arviz     : 0.17.0
preliz    : 0.3.7
matplotlib: 3.8.2

Watermark: 2.4.3