# Lezione 14/04/2026 
# Variabile aleatorie, distribuzioni di probabilità,
# e confronto tra distribuzioni 

# V.A Discrete ----
## V.A. Binomiale ----

# X: v.a. che descrive il numero di Teste 
#    risultanti da venti lanci di un moneta (equilibrata). 

n <- 20 
p <- 0.5

# Ci chiediamo la probabilità di ottenere esattamente 5 teste 

# A mano: choose() permette di ottenere il coefficiente binomiale
choose(n,5) * p^5*(1-p)^(n-5)

# costruendo una funzione 
binomiale <- function(x, n, p){
  choose(n,x) * p^x*(1-p)^(n-x)
}
binomiale(5, n = 20, p = 0.5)

# utilizzando la funzione R
dbinom(5, size = 20, prob = 0.5)

# Esercizio: Calcola la probabilità di ottenere 
# esattamente due volte il 6 lanciando 15 volte 
# un dado regolare.
dbinom(2, size = 15, prob = 1/6) 

# Funzione di probabilità
x <- 0:20
pbin <- binomiale(x,n,p)
sum(pbin)  # Verifica somma a 1 

# Grafico funzione di probabilità
plot(0:20, pbin, type = "h",
     xlab = "x", ylab = "P(X = x)",
     main = "Binomiale(20, 0.5)")

# Funzione di ripartizione 
# Calcoliamo la probabilità di ottenere al più 5 teste in 20 
# lanci di una moneta equilibrata
pbinom(5, size = 20, prob = 0.5)
sum(dbinom(0:5, size = 20, prob = 0.5))

# e la probabilità di ottenere più di 5 teste in 20 
# lanci di una moneta equilibrata
pbinom(5, size = 20, prob =  0.5, lower.tail = F)
1-pbinom(5, size = 20, 0.5)

# Grafico Funzione di Ripartizione 
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: Calcola la probabilità di ottenere 
# almeno 7 volte il 6 lanciando 15 volte un dado regolare.
pbinom(6, size = 15, prob = 1/6) 

## V.A. Poisson ----
# 1. Funzione di probabilità Poisson
poisson_pmf <- function(x, lambda) {
  exp(-lambda) * lambda^x / factorial(x)
}

# Confronto con dpois()
poisson_pmf(3, 5)
dpois(3, 5)

# 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))

# V.A Continue ----
## V.A. Normale ----

# Funzione di Densità 
dnorm(1,mean = 1, sd = 2)
xx <- seq(-6,6, 0.01)
# Tre Normali: N(-1,1); N(0,1); N(2,1)
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))

# Tre normali N(0, 1/9); N(0,1); N(0, 3)
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))


# Funzione di ripartizione:
pnorm(1,mean = 1, sd = 2)

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

# Tre Normali: N(-2,1); N(0,1); N(0,9)
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))


# Alternativa a plot()+points(): curve()
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))


# Esercizio: X ~ N(170, sigma^2 = 100) 
# Si calcoli la P(X <= 185) e la P(165 <= X <= 190)
pnorm(185, mean = 170, sd = 10)
pnorm(190, mean = 170, sd = 10) - pnorm(165, mean = 170, sd = 10)

## V.A. Uniforme ----
x <- seq(0,1,0.1)
dunif(x)  #valore costante in ogni punto del supporto 
punif(x)  # crescita lineare

# 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: X ~ U(0, 2)
# Si calcoli la P(X <= 1.5) e P(1 <= X <= 1.5)$
punif(1.5, 0, 2)
punif(1.5, 0, 2) - punif(1, 0, 2) 
# Equivalentemente base x altezza
(1.5 - 1) * 0.5 

# Funzione quantile ----
## Normale ----
# Quantile di ordine 0.25 di una normale di media 1 
# e varianza 4 è
qnorm(0.25, mean = 1, sd = 2)
  
# Grafico 
# 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")

## Binomiale ----
# Quantile di ordine 0.25 di una binomiale con n = 20 e p=0.5 
qbinom(0.25, size = 20, prob= 0.5)

# Grafico
p <- seq(0.001, 0.999, length=1000)
q <- qbinom(p, size = 20, prob = 0.5) 

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

# Simulazione di numeri casuali ----
## Binomiale ----
# Simuliamo 100 lanci da una moneta equilibrata
lanci <- rbinom(100, size = 1, prob = 0.5)
lanci
table(lanci)              # frequenze osservate
prop.table(table(lanci))   # proporzioni

## Poisson ----
# simuliamo 500 numeri da una poisson con lambda = 3
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))

# Uniforme ----
# Simuliamo 55 numeri a una uniforme(0,1)
set.seed(14)
simu <- runif(500, min = 0, max = 1)

sum(simu < 0.4) / length(simu)  #proporzione osservata 
punif(0.4, 0, 1)                # probabilità teorica

# Normale ----
# Simuliamo 1000 numeri da una N(10, 25)
set.seed(42)
dati.norm <- rnorm(1000, mean = 10, sd = sqrt(25))

sum(dati.norm <= 5) / length(dati.norm)
pnorm(5, mean = 10, sd =5)

# 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 ~ Binom(20, 0.3) 
# (media teorica np = 6):
# Quindi:
# 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 <- 20
p <- 0.3

n_vals <- c(10, 50, 100, 500, 1000, 5000)
medie_campionarie  <- sapply(n_vals, 
                             function(x) mean(rbinom(x, size = n, prob = p)))
media_teorica <- n*p
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 = media_teorica, 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))


# Confronto tra dstribuzioni empiriche e teoriche ----
## Confronto tra ecdf ----
# Genero da Normale
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)

# Genero da t-student(df = 2): distribuzione 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)

## Confronto tra Densità ----
# Istogramma verso funzione di densità
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))

# Funzioni densità kernel 
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))

## Grafico Quantile-Quantile 
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))

# Se i dati non vengono da una 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)

# Interpolazione: confronto tra numerosità diverse
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)")


# Confronto tra due insiemi di dati osservati 
data(iris)
virgi  <- iris$Sepal.Length[iris$Species == "virginica"]
versi  <- iris$Sepal.Length[iris$Species == "versicolor"]
setosa <- iris$Sepal.Length[iris$Species == "setosa"]

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

## Confronto delle Funzioni di Ripartizione Empiriche
  
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)

## Confronto Istogrammi e Funzioni di Densità
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)

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))

# In finestra grafica 2 x 2
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 Pacchetto `lattice`
library(lattice)
histogram(~ Sepal.Length | Species, data = iris,
          type = "density", layout = c(1, 3))
densityplot(~ Sepal.Length | Species, data = iris, layout = c(1, 3))

## Il Grafico Quantile-Quantile fra Due Campioni
length(setosa); length(versi)   # entrambi n = 50

# A mano
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
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)

# Usando qqplot
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