Funzioni legate alle variabili casuali in R

Le quattro funzioni per le variabili casuali

In R sono disponibili una serie di funzioni per lavorare con i principali modelli probabilistici ed altre che seguono la stessa logica di funzionamento sono disponibili in alcuni package aggiuntivi. Le funzioni sono tutte caratterizzate da uno stesso prefisso (la prima lettere del nome) che ne identifica l’obiettivo, e da un suffisso che richiama il particolare modello probabilistico:

• r_ per la generazione di numeri (pseudo) casuali dal particolare modello

• d_ per la distribuzione di probabilità nel caso di modelli discreti / per la funzione di densità di probabilità nel caso di modelli continui

• p per la funzione di ripartizione

• q per la funzione quantile, ovvero per l’inversa della funzione quantile.

Nei due paragrafi seguenti mostro un esempio per un modello discreto (la v.c. binomiale) e per un modello continuo (la v.c. normale).

La variabile casuale binomiale

Il suffisso per lavorare con un modello binomiale è binom, per cui le funzioni di interesse sono:

Funzione	Obiettivo
dbinom()	Distribuzione di probabilità
pbinom()	Funzione di ripartizione
qbinom()	Funzione quantile
rbinom()	Generazione di numeri casuali

La funzione dbinom

Il seguente codice calcola la probabilità per un dato valore di interesse (numero di successi osservati in n prove):

dbinom(x = 5, prob = 1/6, size = 10)

[1] 0.01302381

Facendo riferimento ad un esempio concreto, si tratta della probabilità di osservare 5 volte una faccia di interesse (il successo, che ha appunto probabilità 1/6) in 10 lanci del dado. Usando la logica vettoriale di R è facile calcolare la distribuzione di probabilità completa della variabile in questione passando in input l’intero supporto della variabile:

dbinom(x = 0:10, prob = 1/6, size = 10)

 [1] 1.615056e-01 3.230112e-01 2.907100e-01 1.550454e-01 5.426588e-02
 [6] 1.302381e-02 2.170635e-03 2.480726e-04 1.860544e-05 8.269086e-07
[11] 1.653817e-08

Per rendere i numeri più leggibili è possibile evitare la notazione esponenziale di default chiedendo un prefissato numero di cifre decimali:

round(dbinom(x = 0:10, prob = 1/6, size = 10), 6)

 [1] 0.161506 0.323011 0.290710 0.155045 0.054266 0.013024 0.002171 0.000248
 [9] 0.000019 0.000001 0.000000

La verifica empirica che si tratta di una distribuzione di probabilità è presto fatta:

sum(dbinom(x = 0:10, prob = 1/6, size = 10))

[1] 1

Possiamo ottenere una rappresentazione grafica del modello in questione usando la grafica di base (funzione plot, in attesa di mostrare l’altro sistema grafico che utilizziamo per le analisi su dati reali:

ascisse <- 0:10
ordinate <- dbinom(x = ascisse, prob = 1/6, size = 10)
plot(ascisse, ordinate, type = "h")

E’ chiaramente possibile evitare il passaggio per i due oggetti ascisse ed ordinate scrivendo gli input opportuni direttamente nella funzione plot

plot(0:10, dbinom(x = 0:10, prob = 1/6, size = 10), type = "h")

Allo stesso modo è possibile provare a “parametrizzare” il codice in modo che sia veloce esplorare il comportamento della distribuzione binomiale al variare dei parametri:

n <- 10
pi_greco <- 1/6
ascisse <- 0:n
ordinate <- dbinom(x = ascisse, prob = pi_greco, size = n)
plot(ascisse, ordinate, type = "h")

Questo consente di cambiare solo i valro di n e pi_greco per vedere come cambia la distribuzione. In particolare cambiando il valore di pi_greco mantenendo n fisso è possibile studiare la forma della distribuzione al cambiare della probabilità del successo nella singola prova, mentre facendo crescere n mantenendo fisso pi_greco il comportamento al limite (teorema del limite centrale nella versione di De Moivre - Laplace).

La funzione pbinom

La funzione pbinom permette di ottenere la funzione di ripartizione ed impostando ad s l’argomento type si può ottenere facilmente la corrispondente rappresentazione grafica:

plot(0:10, pbinom(q = 0:10, prob = 1/6, size = 10), type = "s")

La funzione qbinom

La funzione quantile è l’inversa della funzione quantile, che restituisce cioè per una data probabilità di interesse (asse orizzontale) il corrispondente quantile della variabile (asse verticale). Ecco una veloce e grossolana rappresentazione per avere un’idea dell’andamento sempre usando lo stesso esperimento binomiale:

plot(seq(0, 1, 0.01), qbinom(p = seq(0, 1, 0.01), prob = 1/6, size = 10), type = "l")

La funzione rbinom`

Infine la funzione rbinom permette di ottenere numeri (pseudo)casuali a partire dal modello considerato. In altre parole la funzione permette di simulare un esperimento binomiale e contare il numero di successi ottenuto:

rbinom(n = 1, size = 10, prob = 1/6)

[1] 1

Il primo argomento n permette di specificare il numero di esperimenti binomiali da generare (numero di osservazioni casuali da generare). Se si vuole simulare per n volte il lancio di 10 dadi (o 10 volte il lancio dello stesso dato) contando il numero di successi in ciascuna prova è sufficiente cambiare il valore dell’argomento di n:

rbinom(n = 100, size = 10, prob = 1/6)

  [1] 4 3 1 2 2 3 2 2 2 0 2 2 2 1 5 1 2 0 3 2 3 2 1 0 3 0 2 0 2 4 2 2 1 3 1 1 0
 [38] 1 1 3 1 1 2 0 1 0 4 2 2 2 1 1 1 2 0 1 2 1 2 1 1 2 1 0 1 1 1 4 1 3 1 1 1 4
 [75] 0 1 0 1 1 1 1 1 3 3 3 1 2 3 0 1 2 1 0 1 3 2 1 1 2 0

Una semplice verifica empirica di quello che sappiamo dalla teoria:

mean(rbinom(n = 1000, size = 10, prob = 1/6))

[1] 1.612

ovvero che il valore atteso della binomiale è \(n \times \pi = 10 \times 0.1666667 = 1.6666667\), che è abbastanza vicino alla media delle 1000 ripetizioni sopra effettuate dell’esperimento binomiale.

Nota

Le funzioni per generare numeri casuali in R (così come in qualunque altro linguaggio) utilizzano un algoritmo deterministico per ottenere la sequenza di numeri. Questi algoritmi partono fissando un numero inziale, detto seme, da cui si ottiene la sequenza di numeri, che vengono pertanto detti pseudo-casuali. Per ottenere la riproducibilità dei risultati (stessi valori dei numeri ottenuti in qualunque esecuzione delllo stesso codice, anche su PC differenti) è possibile fissare il seme dell’algoritmo con la funzione set.seed():

set.seed(seed = 123)
rbinom(5, size = 10, prob = 1/6)

[1] 1 3 1 3 4

set.seed(seed = 123)
rbinom(5, size = 10, prob = 1/6)

[1] 1 3 1 3 4

La strategia comunemente utilizzata per le simulazioni con i numeri casuali consiste nel fissare un seme prima dell’utilizzo della prima chiamata di una funzione che restituisce numeri casuali: tutte le sequenze successive saranno regolate dello stesso seme iniziale:

# fisso un seme iniziale
set.seed(seed = 123)
# genero 3 numeri casuali da un modello binomiale
rbinom(3, size = 10, prob = 1/6)

[1] 1 3 1

# faccio altre operazioni
x <- seq(1:10)
y <- mean(x)
z <- sd(x)
# genero altri 3 numeri casuali da un modello binomiale
rbinom(3, size = 10, prob = 1/6)

[1] 3 4 0

# fisso lo stesso seme iniziale
set.seed(seed = 123)
# genero una prima sequenza di 3 numeri casuali
rbinom(3, size = 10, prob = 1/6)

[1] 1 3 1

# ed una seconda sequenza di 3 numeri casuali
rbinom(3, size = 10, prob = 1/6)

[1] 3 4 0

Dall’esempio si nota come le due sequenze ottenute sono identiche anche se nel primo caso tra le due ho effettuato altri calcoli.

La variabile casuale normale

Il suffisso per lavorare con un modello normale è norm, per cui le funzioni di interesse sono:

Funzione	Obiettivo
dnorm()	Funzione di densità di probabilità
pnorm()	Funzione di ripartizione
qnorm()	Funzione quantile
rnorm()	Generazione di numeri casuali

La funzione dnorm

Il seguente codice calcola il valore della funzione di densità di probabilità (assse verticale) di un modello normale con media \(\mu = 10\) e scarto quadratico medio \(\sigma = 2\) calcolata nel punto \(x = 9\) (asse orizzontale):

dnorm(x = 9, mean = 10, sd = 2)

[1] 0.1760327

Nota

La parametrizzazione “canonica” del modello normale è effettuata usando la varianza \(\sigma^2\) ma la funzione dnorm lavora usando lo scarto quadratico medio \(\sigma\).

Usando una griglia di valori è facile ottenere anche una prima rappresentazione grafica usando sempre la grafica di base di R:

ascisse <- seq(2, 18, 0.01)
ordinate <- dnorm(x = ascisse, mean = 10, sd = 2)
plot(ascisse, ordinate, type = "l", xlab = "x", ylab = "f(x)")

La funzione pnorm

Usando il vettore di ascisse sopra generato è facile calcolare e rappresentare la funzione di ripartizione dello stesso modello:

ordinate_fr <- pnorm(q = ascisse, mean = 10, sd = 2)
plot(ascisse, ordinate_fr, type = "l", xlab = "x", ylab = "F(x)")

La funzione qnorm

Per tracciare la funzione quantile è necessario invece ragionare sulla scala della probabilità ([0, 1]), poichè la funzione quantile restituisce il quantile di un modello in corrispondenza di un dato valore di probabilità:

ascisse_q <- seq(0, 1, 0.01)
ordinate_q <- qnorm(p = ascisse_q, mean = 10, sd = 2)
plot(ascisse_q, ordinate_q, type = "l", xlab = "p", ylab = "Q(p)")

La funzione rnorm

La sintassi della funzione rnormè simile a quella vista per le funzioni precedenti ma cambia per il primo argomento in input che indica il numero di valori (pseudo)casuali che si è interessati a generare:

rnorm(n = 5, mean = 10, sd = 2)

[1] 10.141017 10.258575 13.430130 10.921832  7.469878

Anche in questo caso, come già visto per la funzione dbinom, è conveniente fissare il seme per regolare la sequenza di valori ottenuti:

set.seed(123)
mille_numeri <- rnorm(n = 1000, mean = 10, sd = 2)

Possiamo fare semplici verifiche empiriche per esplorare la sequenza ottenuta:

# un istogramma dei 1000 numeri ottenuti
hist(mille_numeri, main = '1000 "normali" numeri (pseudo)casuali')

# media e scarto quadratico medio
c(media = mean(mille_numeri), sqm = sd(mille_numeri))

   media      sqm 
10.03226  1.98339 

# sintesi a 5 (minimo, massimo e tre quartili)
fivenum(mille_numeri)

[1]  4.380451  8.742515 10.018419 11.329576 16.482080

Altre verifiche, ai fini di esercizi di riepilogo su puntamenti ed operatori logici, possono consistere nel controllare la percentuale di punti in alcuni intervalli prefissati:

# sappiamo (anche consultando le tavole) la percentuale di valori in un intervello centrato sulla media ampio 2 scarti quadratici medi
round(100 * (pnorm(12, mean = 10, sd = 2) - pnorm(8, mean = 10, sd = 2)), 4)

[1] 68.2689

# 4 scarti quadratici medi (2 a sinistra della media e 2 a destra)
round(100 * (pnorm(14, mean = 10, sd = 2) - pnorm(6, mean = 10, sd = 2)), 4)

[1] 95.45

# 6 scarti quadratici medi (3 a sinistra della media e 3 a destra)
round(100 * (pnorm(16, mean = 10, sd = 2) - pnorm(4, mean = 10, sd = 2)), 4)

[1] 99.73

Possiamo verificarlo empiricamente sui 1000 numeri ottenuti:

# numero di valori che soddisfa la condizione di essere compreso nell'intervallo [8, 12] (condizione in AND, operatore &)
sum(mille_numeri <= 12 & mille_numeri >= 8)

[1] 678

# per calcolare la percentuale basta dividere semplicemente per il numero di valori e moltiplicare per 100
100 * sum(mille_numeri <= 12 & mille_numeri >= 8) / length(mille_numeri)

[1] 67.8

Stesso “controllo” può essere fatto per gli altri due intervalli considerati:

# ecco la percentuale di valori tra -2 scarti e + 2 scarti dalla media
100 * sum(mille_numeri <= 14 & mille_numeri >= 6) / length(mille_numeri)

[1] 95.3

# e tra -3 scarti e +3 scarti dalla media
100 * sum(mille_numeri <= 16 & mille_numeri >= 4) / length(mille_numeri)

[1] 99.9