1 Distribuzioni di Probabilità in R

R dispone di funzioni predefinite per le principali distribuzioni di probabilità, sia discrete che continue. Per ciascuna famiglia il nome segue una convenzione fissa: un prefisso di una lettera seguito dal nome della distribuzione.

Prefisso Funzione calcolata Esempio (norm)
d funzione di probabilità / densità dnorm(x, mean, sd)
p funzione di ripartizione pnorm(x, mean, sd)
q funzione quantile qnorm(p, mean, sd)
r generazione di valori casuali rnorm(n, mean, sd)

R implementa le principali famiglie di distribuzioni di probabilità con la convenzione: prefisso d/p/q/r + nome della distribuzione.

Distribuzione Nome R Parametri aggiuntivi
Normale norm mean, sd
t di Student t df
Chi-quadrato chisq df
F di Fisher f df1, df2
Gamma gamma shape, scale
Uniforme unif min, max
Binomiale binom size, prob
Poisson pois lambda
Esponenziale exp rate
Log-normale lnorm meanlog, sdlog
Weibull weibull shape, scale

💡 Per l’elenco completo: ?Distributions

2 Distribuzioni Discrete

2.1 Distribuzione Binomiale

La distribuzione Binomiale \(X \sim \text{Binom}(n, p)\) modella il numero di successi in \(n\) prove indipendenti, dove per ogni prova la probabilità dell’evento successo è pari a \(p\).

La funzione di probabilità è: \[P(X = x) = \binom{n}{x} p^x (1-p)^{n-x}, \quad x = 0, 1, \ldots, n\]

2.1.1 Calcolo Manuale e con le Funzioni R

Supponiamo che \(X\) sia la variabile aleatoria che descrive il numero di Teste risultanti da venti lanci di un moneta (equilibrata).

In questo caso:

  • \(n = 20\) prove indipendenti (il risultato del secondo lancio non dipende da cosa sia successo al primo, e coì per gli altri)
  • \(p = 0.5\) probabilità di successo (in questo caso assumiamo che sia l’uscita di Testa)
  • \(X\sim Bin(n=20, p = 0.5)\)

Calcoliamo la probabilità che di avere esattamente 5 teste, ovvero \(P(X=5)\)

  • a mano
choose(20, 5) * 0.5^5 * (1 - 0.5)^15
#> [1] 0.01478577
  • a mano con funzione personalizzata
binomiale <- function(x, n, p) {
  choose(n, x) * p^x * (1 - p)^(n - x)
}
binomiale(5, 20, 0.5)
#> [1] 0.01478577
  • usando dbinom()
dbinom(5, size = 20, prob = 0.5)
#> [1] 0.01478577

📌 Esercizio veloce: calcola la probabilità di ottenere esattamente due volte il 6 lanciando 15 volte un dado regolare.

\(X \sim \text{Binom}(15,\, 1/6)\)

binomiale(2, 15, 1/6)
#> [1] 0.272603
dbinom(2, size = 15, prob = 1/6)   # verifica
#> [1] 0.272603

2.1.2 Funzione di Probabilità

Calcoliamo \(P(X = x)\) per tutti i valori \(x = 0, 1, \ldots, 20\) e verifichiamo che la somma sia 1, condizione affinchè \(P(X = x)\) sia una distribuzione di probabilità.

\[\sum_{x=0}^{n}\binom{n}{x} p^x (1-p)^{n-x} = 1\] Utilizziamo la funzione appena costruita

x    <- 0:20
pbin <- binomiale(x, 20, 0.5)
pbin
#>  [1] 9.536743e-07 1.907349e-05 1.811981e-04 1.087189e-03 4.620552e-03
#>  [6] 1.478577e-02 3.696442e-02 7.392883e-02 1.201344e-01 1.601791e-01
#> [11] 1.761971e-01 1.601791e-01 1.201344e-01 7.392883e-02 3.696442e-02
#> [16] 1.478577e-02 4.620552e-03 1.087189e-03 1.811981e-04 1.907349e-05
#> [21] 9.536743e-07
sum(pbin)   # deve essere 1
#> [1] 1

Visualizziamo il grafico della funzione di probabilità

plot(0:20, pbin, type = "h",
     xlab = "x", ylab = "P(X = x)",
     main = "Binomiale(20, 0.5)")

💡 Utilizziamo type = "h" che ci permette di diegnare dei bastoncini.

2.1.3 Funzione di Ripartizione

La funzione di ripartizione è per definizione \[F_X(x) = P(X \leq x)\] Nel caso della binomiale si ottiene come \[F_X(x)= \sum_{k=0}^{x}P(X=k) = \sum_{k=0}^{x} \binom{n}{k} p^k (1-p)^{n-k}\]

In R otteniamo i valori della funzione di ripartizione calcolati in x mediante pbinom.

Quindi se si vuole ottenere la probabilità:

  • di osservare al più 5 successi (in \(n\) lanci con probabilità di successo pari a \(p\)): \(P(X\leq 5)\)
pbinom(5, 20, 0.5) 
#> [1] 0.02069473
sum(binomiale(0:5,20,0.5))
#> [1] 0.02069473
  • di osservare più di 5 successi: \(P(X > 5) = 1 - P(X\leq 5)\)
1-pbinom(5, 20, 0.5) 
#> [1] 0.9793053
# Equivalentemente
pbinom(5, 20, 0.5, lower.tail = F)
#> [1] 0.9793053

Trattandosi di una distribuzione discreta il grafico della funzione di ripartizione (in inglese , cdf) è una funzione a gradini.

cdf <- pbinom(0:20, 20, 0.5) 
plot(0:20, cdf)
segments(x[-length(x)], cdf[-length(cdf)],
         x[-1], cdf[-length(cdf)])
segments(x, c(0, cdf[-length(cdf)]),
         x, cdf, lty = "dashed")

📌 Esercizio veloce: calcola la probabilità di ottenere almeno 7 volte il 6 lanciando 15 volte un dado regolare.

\(X \sim \text{Binom}(15,\, 1/6)\)

1-sum(binomiale(0:6, 15, 1/6))
#> [1] 0.006603585
1 - pbinom(6, 15, 1/6)   # 
#> [1] 0.006603585
pbinom(6, 15, 1/6, lower.tail = F)
#> [1] 0.006603585

2.1.4 Esercizio

Esercizio 1 — Distribuzione di Poisson La distribuzione di Poisson \(X \sim \text{Pois}(\lambda)\) modella il numero di eventi rari in un intervallo fissato. La funzione di probabilità è: \[P(X = x) = \frac{e^{-\lambda} \lambda^x}{x!}, \quad x = 0, 1, 2, \ldots\]

  1. Scrivi una funzione poisson_pmf(x, lambda) che valuti la funzione di probabilità
  2. Rappresenta graficamente la funzione di probabilità per \(\lambda = 3\)
  3. Extra (per casa): dimostra graficamente che una \(\text{Binom}(n, p)\) con \(n\) grande e \(p\) piccolo è ben approssimata da una \(\text{Pois}(\lambda = np)\). Usa ad esempio \(n = 100\), \(p = 0.02\).
# 1. Funzione di probabilità Poisson
poisson_pmf <- function(x, lambda) {
  exp(-lambda) * lambda^x / factorial(x)
}

# Confronto con dpois()
poisson_pmf(3, 5)
#> [1] 0.1403739
dpois(3, 5)
#> [1] 0.1403739
# 2. Grafico
x <- 0:15
plot(x, dpois(x, lambda = 3), type = "h",
     xlab = "x", ylab = "P(X = x)",
     main = "Poisson(lambda = 3)")

# 3. Approssimazione Binomiale → Poisson
n <- 100; p <- 0.02; lambda <- n * p
x <- 0:10
plot(x, dbinom(x, n, p), type = "h",
     xlab = "x", ylab = "Probabilità",
     main = paste0("Binom(", n, ",", p, ") vs Pois(", lambda, ")"))
points(x + 0.1, dpois(x, lambda), type = "h", col = 2, lwd = 3)
legend("topright", legend = c("Binomiale", "Poisson"),
       col = c(1, 2), lwd = c(1, 3))

3 Distribuzioni Continue

3.1 Distribuzione Normale

La distribuzione Normale \(X \sim N(\mu, \sigma^2)\) è la distribuzione continua più importante in statistica. È completamente determinata da media \(\mu\) e varianza \(\sigma^2\) (o deviazione standard \(\sigma\)).

⚠️ In R le funzioni dnorm(), pnorm(), qnorm(), rnorm() usano come secondo e terzo argomento media e deviazione standard (non la varianza). Se non si specificano gli argomenti mean e sd si utilizzerà la normale standard, ovvero una normale di media 0 e varianza 1.

3.1.1 Funzione di Densità e di Ripartizione

La funzione di densità di \(X\sim \mathcal{N}(\mu, \sigma^2)\) è data da \[ f_X(x) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-\frac{(x-\mu)^2}{2\sigma^2}}\] Per una generica v.a. continua (non solo di una v.a. normale) la funzione di densità ha le seguenti proprietà:

  • \(f_X(x) \geq 0\)
  • \(\int_{-\infty}^{\infty} f_X(x) = 1\)

La distribuzione normale standard, solitamente indicata con \(Z \sim \mathcal{N}(0,1)\) si ottiene per standardizzazione, ovvero \(Z = \frac{X-\mu}{\sigma}\), e ha densità

\[ f_Z(z) = \frac{1}{\sqrt{2 \pi}} e^{-\frac{z^2}{2}}\]

In R, assumendo una normale di media 1 e varianza 4, nel punto 1 vale

# Calcoliamo nel punto 1. 
dnorm(1, mean = 1, sd = 2)
#> [1] 0.1994711

Quindi possiamo rappresentarla:

  • al variare di \(\mu\) (fissato \(\sigma\)) si ha
xx <- seq(-6, 6, by = 0.01)

plot(xx, dnorm(xx, mean = -1, sd = 1), type = "l",
     xlim = c(-6, 6), ylim = c(0, 0.45),
     xlab = "x", ylab = "f(x)",
     main = "Densità Normale al variare di mu")
points(xx, dnorm(xx, mean = 0,  sd = 1), 
       type = "l", lty = 2, col = 2, lwd = 3)
points(xx, dnorm(xx, mean = 2,  sd = 1), 
       type = "l", lty = 3, col = 3, lwd = 3)
legend("topright",
       legend = c("N(-1, 1)", "N(0, 1)", "N(2, 1)"),
       col = c(1, 2, 3), lty = 1:3, lwd = c(1, 3, 3))

  • Al variare di \(\sigma\) (fissato \(\mu\)) si ha
xx <- seq(-6, 6, by = 0.01)

plot(xx, dnorm(xx, mean = 0, sd = 1/3), type = "l",
     xlim = c(-6, 6), ylim = c(0, 1.45),
     xlab = "x", ylab = "f(x)",
     main = "Densità Normale al variare di sigma")
points(xx, dnorm(xx, mean = 0,  sd = 1), 
       type = "l", lty = 2, col = 2, lwd = 3)
points(xx, dnorm(xx, mean = 0,  sd = 3), 
       type = "l", lty = 3, col = 3, lwd = 3)
legend("topright",
       legend = c("N(0, 1/9)", "N(0, 1)", "N(0, 9)"),
       col = c(1, 2, 3), lty = 1:3, lwd = c(1, 3, 3))

La funzione di ripartizione per una generica v.a. è data da \[F_X(x) = P(X\leq x) = \int_{-\infty}^{x}f_X(x)dx\] con

  • \(0 \leq F_x(x) \leq 1\)
  • \(P(a \leq X \leq b) = P(X \leq a) - P(X \leq b) = F_X(a) - F_X(b)\)

In R, possiamo ottenere il valore della funzione di ripartizione in un punto (previa specificazione di media e deviazione standard). Ad esempio nel punto 1, per una normal edi media 1 e varianza 4, si ha

# Calcoliamo nel punto 1. 
pnorm(1, mean = 1, sd = 2)
#> [1] 0.5

Quindi possiamo rappresentarla per diversi valori di \(\mu\) e \(\sigma\).

xx <- seq(-6, 6, by = 0.01)

plot(xx, pnorm(xx, mean = -2, sd = 1), type = "l",
     xlim = c(-6, 6), 
     xlab = "x", ylab = "f(x)",
     main = "Funzione di Ripartizione")
points(xx, pnorm(xx, mean = 0,  sd = 1), 
       type = "l", lty = 2, col = 2, lwd = 3)
points(xx, pnorm(xx, mean = 0,  sd = 3), 
       type = "l", lty = 3, col = 3, lwd = 3)
legend("topright",
       legend = c("N(-2, 1)", "N(0, 1)", "N(0, 3)"),
       col = c(1, 2, 3), lty = 1:3, lwd = c(1, 3, 3))

💡 La funzione curve() è un’alternativa a plot() + points() per tracciare funzioni matematiche: accetta direttamente un’espressione in x senza dover definire una griglia di punti.

par(mfrow=c(1,2))
# Densità con curve()
curve(dnorm(x, mean = 0, sd = 1), col = 2, lwd = 3, lty = 2,
      xlim = c(-5, 5), xlab = "x", ylab = "f(x)", ylim = c(0,0.6),
      main = "Funzione di Densità")
curve(dnorm(x, mean = -1, sd = 2), col = 4, lwd = 3, lty = 4,
      add = TRUE)
legend("topright",
       legend = c("N(0, 1)", "N(-1, 2)"),
       col = c(2, 4), lty = c(2, 4), lwd = c(3, 3))

# Funzione di ripartizione
curve(pnorm(x, mean = 0, sd = 1), col = 2, lwd = 3, lty = 2,
      xlim = c(-4, 4), xlab = "x", ylab = "F(x)",
      main = "Funzione di ripartizione")
curve(pnorm(x, mean = -1, sd = 2), col = 4, lwd = 3, lty = 4,
      add = TRUE)
legend("topleft",
       legend = c("N(0, 1)", "N(-1, 4)"),
       col = c(2, 4), lty = c(2, 4), lwd = c(3, 3))

par(mfrow=c(1,1))

3.1.2 Esercizio

Esercizio veloce — Calcoli con la Normale
Sia \(X \sim N(170, 100)\) (media 170, varianza 100, quindi \(\sigma = 10\)).

Calcola:

  1. \(P(X \leq 185)\)
  2. \(P(165 \leq X \leq 190)\)
mu <- 170; sigma <- sqrt(100)

# 1. P(X <= 185)
pnorm(185, mu, sigma)
#> [1] 0.9331928
# 2. P(165 <= X <= 190)
pnorm(190, mu, sigma) - pnorm(165, mu, sigma)
#> [1] 0.6687123

3.2 Distribuzione Uniforme

La distribuzione Uniforme \(X \sim U(a, b)\) è completamente determinata dal minimo valore del supporto (\(a\)) e dal massimo (\(b\)). La sua densità è \[f_X(x) = \frac{1}{b-a} \quad \text{se} \quad a \leq x \leq b \] e la sua funzione di ripartizione è \[F_X(x) = \begin{cases} 0 \quad \text{se}\quad x \leq a \\ \frac{x-a}{b-a} \quad \text{se}\quad a < x \leq b \\ 1 \quad \text{se}\quad x > b \end{cases} \]

⚠️ In R le funzioni dunif(), punif(), qunif(), runif(). Se non si specificano gli argomenti min e max si utilizzerà la uniforme (standard), ovvero una uniforme in (0,1).

Per una distribuzione uniforme continua \(X \sim U(a,b)\) vale che se \(b_1 - a_1 = c_1\) e \(b_2-a_2=c_2\) con \(c_1 = c_2 = c\) allora
\[P(a_1 \leq X \leq b_1) = P(a_2 \leq X \leq b_2) = c \times \frac{1}{b-a}\]

x <- seq(0,1,0.1)
dunif(x)  #valore costante in ogni punto del supporto 
#>  [1] 1 1 1 1 1 1 1 1 1 1 1
punif(x)  # crescita lineare
#>  [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
# Considerando l'uniforme in (0,2)
a <- 0; b <- 2

par(mfrow=c(1,2))
curve(dunif(x, min = a, max = b), 
      xlim = c(-0.2,2.2), col = 4, lwd = 2,
      main = "Densità U(0,2)")

curve(punif(x, min = a, max = b), 
      xlim = c(-0.2,2.2), col = 4, lwd = 2,
      main = "Ripartizione U(0,2)")

par(mfrow=c(1,1))

Esercizio veloce — Calcoli con l’uniforme
Sia \(X \sim U(0, 2)\).

Calcola:

  1. \(P(X \leq 1.5)\)
  2. \(P(1 \leq X \leq 1.5)\)
# Considerando l'uniforme in (0,1)
a <- 0; b <- 2

# 1. 
punif(1.5, min = 0, max = 2)
#> [1] 0.75
#2.
punif(1.5, min = 0, max = 2) - punif(1, min = 0, max = 2)
#> [1] 0.25
# Equivalentemente base x altezza
(1.5 - 1) * 0.5 
#> [1] 0.25

4 La funzione quantile

Per un dato valore \(p \in (0,1)\) la funzione quantile \(Q(p)\) restituisce il valore \(x_p\) appartenente al supporto della variabile aleatoria \(X\). Per v.a.:

  • Discrete: \(x_p = Q(p) = \min{x: F(x) \geq p}\) (non è una inversa pura)
  • Continue: \(x_p = Q(p) = F^{-1}(p)\)

4.1 Normale

👉 qnorm(p) restituisce il quantile della normale (Standard), cioè il valore sotto cui cade la proporzione \(p\) dei dati.

Ad esempio il quantile di ordine \(0.25\) di una normale di media 1 e varianza 4 è:

qnorm(0.25, mean = 1, sd = 4)
#> [1] -1.697959

Possiamo rappresentare la funzione quantile in questo modo

# Considerando la normale standard
p <- seq(0.001, 0.999, length=1000)
q <- qnorm(p)

plot(p, q, type="l",
     xlab="p",
     ylab="Q(p)",
     main="Funzione quantile della Normale standard")

4.2 Binomiale

Nel caso di una v.a. discreta la funzione quantile non sarà continua ma sarà a gradini.

Ad esempio il quantile di ordine 0.25 di \(X \sim Bin(n=20, p = 0.5)\) è dato da

qbinom(0.25, size = 20, prob= 0.5)
#> [1] 8
p <- seq(0, 1, length=1000)
q <- qbinom(p, size = 20, prob = 0.5) 

plot(p, q, type="l",
     xlab="p",
     ylab="Q(p)",
     main="Funzione quantile della Binomiale (n = 20, p = 0.5)")

5 Generazione di Numeri Casuali

R genera numeri pseudo-casuali: sequenze deterministiche che si comportano statisticamente come se fossero casuali. La funzione set.seed() fissa il punto di partenza della sequenza, garantendo la riproducibilità dei risultati.

💡 Usare sempre set.seed() prima di qualsiasi simulazione che si vuole riprodurre esattamente.

5.1 Simulazione di Lanci di Moneta

rbinom() genera valori casuali dalla distribuzione Binomiale. Con \(n = 1\) (size = 1) si ottiene una Bernoulli.

Consideriamo 100 lanci di una moneta equilibrata (0 = croce, 1 = testa), pertanto specifichiamo \(p=0.5\) (prob=0.5).

lanci <- rbinom(100, size = 1, prob = 0.5)
lanci
#>   [1] 0 1 1 0 1 0 0 1 0 1 1 0 1 1 1 0 1 0 0 1 1 0 1 1 1 1 1 0 0 0 0 0 1 0 1 1 0
#>  [38] 1 1 0 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 1 0 0 1 1 0 1 0 0 0 0 1 0 0 0 1 0 0 0
#>  [75] 1 1 1 1 1 1 1 0 0 1 1 0 0 0 0 0 0 1 0 0 0 1 1 0 0 0
table(lanci)              # frequenze osservate
#> lanci
#>  0  1 
#> 57 43
prop.table(table(lanci))   # proporzioni
#> lanci
#>    0    1 
#> 0.57 0.43

5.2 Distribuzione di Poisson

Generiamo i valori da una Poisson con parametro \(\lambda\) pari a 3, e confrontiamo le frequenze relative osservate con le probabilità teoriche.

set.seed(3)
n         <- 500
dati.pois <- rpois(n, lambda = 3)

tab  <- table(dati.pois)          # distribuzione di frequenze assolute
tabo <- as.data.frame(tab)

punti <- as.numeric(levels(tabo$dati.pois))  # valori osservati
freq  <- as.numeric(tab / n)                 # frequenze relative

# Grafico: frequenze osservate (nero) vs probabilità teoriche (rosso)
plot(punti, freq, type = "h",
     ylab = "Probabilità", xlab = "x",
     xlim = c(0, 10),
     main = "Poisson(5): simulato vs teorico")
points(punti + 0.1, dpois(punti, lambda = 3),
       type = "h", col = 2, lwd = 3)
legend("topright",
       legend = c("Frequenze simulate", "Probabilità teoriche"),
       col = c(1, 2), lwd = c(1, 3))

5.3 Distribuzione Uniforme

Generiamo 500 valori da U(0,1) e verifichiamo la proporzione di valori inferiori a 0.4. Ci aspettiamo circa il 40% dei valori sotto 0.4 (probabilità teorica)

set.seed(14)
simu <- runif(500, min = 0, max = 1)

sum(simu < 0.4) / length(simu)   
#> [1] 0.396
punif(0.4, 0, 1) 
#> [1] 0.4

5.4 Distribuzione Normale

Generiamo 1000 valori da N(10, 25) e valutiamo la proporzione di valori inferiori a 5, confrontandola con quella teorica.

set.seed(42)
dati.norm <- rnorm(1000, mean = 10, sd = sqrt(25))

sum(dati.norm <= 5) / length(dati.norm)
#> [1] 0.174
pnorm(5, mean = 10, sd =5)
#> [1] 0.1586553

5.4.1 Esercizio

Esercizio 2 — Legge dei Grandi Numeri
La legge dei grandi numeri afferma che la media campionaria converge alla media teorica all’aumentare di \(n\).

Verifica empiricamente questo risultato per \(X \sim \text{Binom}(20, 0.3)\) (media teorica \(= np = 6\)):

  1. Genera campioni di dimensione crescente: \(n = 10, 50, 100, 500, 1000, 5000\)
  2. Calcola la media campionaria per ciascun campione
  3. Rappresenta graficamente come la media campionaria si avvicina al valore teorico
set.seed(123)
n_vals <- c(10, 50, 100, 500, 1000, 5000)
medie_campionarie <- sapply(n_vals, function(n) mean(rbinom(n, size = 20, prob = 0.3)))

plot(n_vals, medie_campionarie, type = "b", log = "x",
     xlab = "n (scala logaritmica)", ylab = "Media campionaria",
     main = "Legge dei Grandi Numeri — Binom(20, 0.3)",
     ylim = c(5, 7))
abline(h = 20 * 0.3, col = 2, lwd = 2, lty = 2)
legend("topright", legend = c("Media campionaria", "Media teorica (6)"),
       col = c(1, 2), lty = c(1, 2), lwd = c(1, 2))

6 Confronto fra Distribuzioni Empiriche e Teoriche

Nell’analisi di variabili quantitative è spesso utile verificare se i dati osservati si conformano a un modello teorico noto — ad esempio la distribuzione Normale. Il problema si può affrontare in modo grafico (approccio esplorativo) o con test statistici formali (inferenza). I due approcci sono tipicamente usati in modo complementare.

6.1 Confronto tra ECDF Empirica e Teorica

Un primo strumento grafico consiste nel sovrapporre la funzione di ripartizione empirica (ecdf()) alla funzione di ripartizione teorica (pnorm(), pt(), ecc.). Se le due curve sono vicine, i dati sono compatibili con quel modello.

set.seed(23)
dat1 <- rnorm(50)   # dati generati da una gaussiana standard

plot(ecdf(dat1), cex = 0.5,
     main = "ECDF empirica vs teorica gaussiana")
curve(pnorm(x), add = TRUE, col = 2, lwd = 2)

Proviamo con dati generati da una t di Student con 2 gradi di libertà — distribuzione simile alla gaussiana ma con code più pesanti:

dat2 <- rt(50, df = 2)

plot(ecdf(dat2), cex = 0.5, verticals = TRUE,
     main = "ECDF da t(2) vs teorica gaussiana")
curve(pnorm(x), add = TRUE, col = 2, lwd = 2)

📌 Nel secondo caso lo scostamento fra le due curve è evidente nelle code, dove la t produce valori più estremi della gaussiana. La distanza massima fra le due curve è la statistica di Kolmogorov-Smirnov: \[D = \sup_x |\hat{F}(x) - F(x)|\]

6.2 Confronto tra Densità Empirica e Teorica

Un secondo approccio consiste nel sovrapporre la densità teorica all’istogramma (con freq = FALSE) o alla curva di densità kernel.

par(mfrow = c(1, 2))

dat3 <- rnorm(100)
hist(dat3, prob = TRUE, xlim = c(-5, 5), ylim = c(0, 0.6),
     main = "n = 100", xlab = "x", col = "lightblue", border = "white")
curve(dnorm(x, 0, 1), col = 2, lwd = 2, add = TRUE)

dat32 <- rnorm(1000)
hist(dat32, prob = TRUE, xlim = c(-5, 5), ylim = c(0, 0.6),
     main = "n = 1000", xlab = "x", col = "lightblue", border = "white")
curve(dnorm(x, 0, 1), col = 2, lwd = 2, add = TRUE)

par(mfrow = c(1, 1))

⚠️ Con \(n = 100\) l’istogramma è instabile e la somiglianza con il modello è solo approssimativa. Con numerosità maggiori la convergenza migliora sensibilmente. Non vanno sopravvalutate le differenze grafiche con pochi dati.

Con la curva di densità kernel il confronto è più fluido:

dat4 <- rt(1000, df = 2)   # t con 2 gdl: code molto più pesanti

plot(density(dat4), xlim = c(-6, 6), ylim = c(0, 0.6),
     main = "Densità kernel da t(1) vs gaussiana")
curve(dnorm(x, 0, 1), col = 2, lwd = 2, add = TRUE)
legend("topright", legend = c("t(2) empirica", "N(0,1) teorica"),
       col = c(1, 2), lwd = c(1, 2))

📌 Le code della t(2) sono molto più pesanti di quelle gaussiane: la curva empirica è più alta ai margini e più bassa al centro.

6.3 Il Grafico Quantile-Quantile (QQ-plot)

Il confronto dei quantili empirici con quelli teorici è lo strumento più potente per valutare l’aderenza a un modello. L’idea: se i dati provengono da una distribuzione \(F\), allora il quantile empirico di ordine \(p\) dovrebbe essere vicino al quantile teorico \(F^{-1}(p)\). Se i due insiemi di quantili si dispongono lungo la retta bisettrice, i dati sono compatibili con il modello.

Usando la convenzione \(p_i = (i - 0.5)/n\):

set.seed(2222)
dat5  <- rnorm(100)
dat5s <- sort(dat5)
qe    <- ((1:100) - 0.5) / 100

par(mfrow = c(1, 2))

# Grafico quantile-quantile costruito manualmente
plot(qnorm(qe), dat5s,
     main = "QQ-plot manuale", xlab = "Quantili teorici N(0,1)",
     ylab = "Quantili empirici")
abline(0, 1, col = 2, lwd = 2)

# Con la funzione qqnorm()
qqnorm(dat5s, main = "qqnorm()")
qqline(dat5s, col = 2, lwd = 2)

par(mfrow = c(1, 1))

💡 qqnorm() è un caso particolare di qqplot() specializzato per il confronto con la gaussiana. qqline() aggiunge la retta di riferimento passante per il primo e terzo quartile.

Cosa succede con dati che non provengono dalla gaussiana?

set.seed(1111)
dat6 <- rt(100, df = 1)
qqnorm(dat6, main = "QQ-plot per dati da t(1)")
qqline(dat6, col = 2, lwd = 2)

📌 I punti alle code si discostano marcatamente dalla retta di riferimento: la t(1) ha code molto più pesanti della gaussiana, e il QQ-plot lo evidenzia con chiarezza.

Il QQ-plot funziona anche con campioni di numerosità diversa: qqplot() interpola i quantili automaticamente.

set.seed(42)
xz <- rnorm(50, mean = 5, sd = 2)
xc <- rchisq(100, df = 5)
qqplot(xz, xc, pch = 19,
       main = "QQ-plot: N(5,4) vs Chi²(5) — distribuzioni diverse",
       xlab = "Quantili N(5,4)", ylab = "Quantili Chi²(5)")

7 Confronto fra Due Insiemi di Dati Osservati

I grafici introdotti per le variabili quantitative e gli strumenti per il confronto fra distribuzioni empirica e teorica si possono usare anche per confrontare due o più insiemi di dati osservati, affiancando o sovrapponendo le rappresentazioni già viste.

Nel seguito consideremo il dataframe iris che riporta le misurazioni in cm della lunghezza e larghezza dei sepali e dei petali di 50 fiori di tre specie diverse di iris.

Consider-iamo, a titolo di esempio, la variabile numerica Sepal.Length ed estraiamo tale misura per ogni specie. Quando effettueremo confronti a coppie consideremo il confronti tra le specie setosa e versicolor.

data(iris)
virgi  <- iris$Sepal.Length[iris$Species == "virginica"]
versi  <- iris$Sepal.Length[iris$Species == "versicolor"]
setosa <- iris$Sepal.Length[iris$Species == "setosa"]

7.1 Boxplot Affiancati

Il boxplot è particolarmente adatto ai confronti: funziona bene anche con campioni piccoli, dove l’istogramma non è affidabile. Si possono affiancare i boxplot di due gruppi passandoli esplicitamente, oppure usando la notazione formula y ~ gruppo per ottenere automaticamente un pannello per ogni livello del fattore.

boxplot(Sepal.Length ~ Species, data = iris,
        main = "Lunghezza sepalo per specie",
        ylab = "Lunghezza (cm)",
        col  = c("lightblue", "lightgreen", "lightyellow"))

💡 Il simbolo ~ definisce una formula in R: la variabile a sinistra viene studiata condizionatamente alla variabile a destra. L’argomento data evita di dover specificare il nome del data frame ad ogni variabile.

7.2 Confronto delle Funzioni di Ripartizione Empiriche

Considerando solo i gruppi setosa e versicolor, sovrapponiamo le ECDF usando add = TRUE nel secondo plot():

plot(ecdf(setosa), verticals = TRUE, xlim = c(4, 9),
     main = "ECDF: setosa vs versicolor",
     xlab = "Lunghezza sepalo")
plot(ecdf(versi), verticals = TRUE, col = 2, add = TRUE)
legend("topleft", legend = c("setosa", "versicolor"),
       col = c(1, 2, 3), lwd = 2)

📌 La ECDF di versicolor \(\hat{F}_v(x)\) è spostata a destra rispetto a quella di setosa \(\hat{F}_s(x)\): \(\hat{F}_v(x) \leq \hat{F}_s(x)\) per ogni \(x\), indicando che i sepali di versicolor sono sistematicamente più lunghi. La distanza massima fra le due curve è la statistica di Kolmogorov-Smirnov: \[D = \sup_x |\hat{F}_s(x) - \hat{F}_v(x)|\]

7.3 Confronto delle Funzioni di Densità

La sovrapposizione di istogrammi è poco efficace quando le distribuzioni si sovrappongono: le barre si coprono rendendo il confronto difficile.

hist(setosa, xlim = c(4, 8.5), freq = FALSE,
     xlab = "Lunghezza sepalo",
     main = "Sovrapposizione di istogrammi (poco efficace)",
     col  = "lightblue", border = "white")
hist(versi, freq = FALSE, border = 2, density = 4.5, add = TRUE)

Più chiaro è il confronto con le curve di densità kernel:

plot(density(setosa), xlim = c(4, 8.5),
     xlab = "Lunghezza sepalo",
     main = "Densità kernel: setosa vs versicolor")
lines(density(versi), col = 2)
legend(x = 7.5, y = 0.6,
       legend = c("setosa", "versicolor"),
       fill   = c(1, 2))

💡 Le curve di densità kernel sono più efficaci degli istogrammi sovrapposti. Tuttavia sono affidabili soprattutto con gruppi di numerosità sufficientemente elevata (indicativamente \(n > 100\)).

Con par(mfrow) si possono affiancare istogramma e densità per ciascun gruppo, ma le scale degli assi rischiano di essere diverse e rendere difficile il confronto visivo:

par(mfrow = c(2, 2))
hist(setosa, freq = FALSE,
     xlab = "lunghezza sepalo", main = "Istogramma Setosa")
plot(density(setosa), main = "Densità Setosa")
hist(versi, freq = FALSE, border = 2, density = 4.5,
     main = "Istogramma Versicolor")
plot(density(virgi), xlim = c(4, 8.5), col = 2,
     xlab = "lunghezza sepalo", main = "Densità Versicolor")

par(mfrow = c(1, 1))

⚠️ Il confronto è difficile perché le scale sull’asse orizzontale non sono le stesse. Per grafici multipli con scala comune è preferibile usare il pacchetto lattice o ggplot2.

7.4 Il Pacchetto lattice

Il pacchetto lattice produce grafici multipli (trellis) con scala comune automatica, usando la notazione formula y | gruppo:

library(lattice)
histogram(~ Sepal.Length | Species, data = iris,
          type = "density", layout = c(1, 3))

densityplot(~ Sepal.Length | Species, data = iris, layout = c(1, 3))

💡 La sintassi y | gruppo in lattice produce automaticamente un pannello per ogni livello del fattore con scala comune, rendendo il confronto visivo molto più affidabile rispetto a par(mfrow).

7.5 Il Grafico Quantile-Quantile fra Due Campioni

Il confronto dei quantili empirici per corrispondenti valori di \(p\) è uno strumento potente per confrontare due distribuzioni. Se le due distribuzioni sono simili, i quantili di pari ordine \(p\) dovrebbero coincidere e i punti del grafico dovrebbero allinearsi lungo la retta bisettrice.

Quando i due campioni hanno la stessa numerosità \(n\), il grafico è semplicemente uno scatterplot dei due vettori ordinati: il \(k\)-esimo valore ordinato di ciascun gruppo rappresenta il quantile di ordine \(k/n\).

length(setosa); length(versi)   # entrambi n = 50
#> [1] 50
#> [1] 50
plot(sort(setosa), sort(versi), pch = 19,
     main = "QQ-plot: setosa vs versicolor",
     xlab = "Quantili setosa", ylab = "Quantili versicolor")
abline(0, 1, col = "gray")   # bisettrice

📌 I punti non sono sulla bisettrice (i sepali di virginica sono più lunghi), ma sembrano disporsi lungo una retta quasi parallela ad essa: le due distribuzioni hanno forma simile ma media diversa.

Una retta di riferimento più robusta passa per i quartili corrispondenti delle due distribuzioni:

qsetosa <- quantile(setosa, probs = c(.25, .75))
qversi <- quantile(versi, probs = c(.25, .75))
s   <- (qversi[2] - qversi[1]) / (qsetosa[2] - qsetosa[1])
int <- qversi[1] - s * qsetosa[1]

plot(sort(setosa), sort(versi), pch = 19,
     main = "QQ-plot con retta sui quartili",
     xlab = "Quantili setosa", ylab = "Quantili versicolor")
abline(int, s, col = 2, lwd = 2)

La funzione qqplot() produce lo stesso risultato direttamente:

qqplot(setosa, versi, pch = 19,
       main  = "QQ-plot con qqplot()",
       xlab  = "Quantili setosa",
       ylab  = "Quantili versicolor")

abline(0, 1, col = "gray")

qsetosa <- quantile(setosa, probs = c(.25, .75))
s2  <- (qversi[2] - qversi[1]) / (qsetosa[2] - qsetosa[1])
int2 <- qversi[1] - s2 * qsetosa[1]
abline(int2, s2, col = 2, lwd = 2)

# i punti sono lungo una retta non parallela alla bisettrice:
# le distribuzioni sono simili ma con diversa media e varianza

7.5.1 Interpretazione del QQ-plot fra due campioni

Pattern osservato Interpretazione
Punti sulla bisettrice Le due distribuzioni sono uguali
Retta parallela alla bisettrice Stessa forma, medie diverse
Retta non parallela Stessa forma, medie e varianze diverse
Nuvola non lineare Le distribuzioni differiscono anche nella forma