# approssimazione quadratica della log-verosimiglianza di Bernoulli -----------
#
# Obiettivo: mostrare come lo sviluppo in serie di Taylor del 2° ordine
# nel punto di massima verosimiglianza (theta_hat) approssima la
# log-verosimiglianza normalizzata al variare di n e della proporzione
# campionaria osservata.


# funzione di log-verosimiglianza normalizzata --------------------------------

# loglik_bern calcola la log-verosimiglianza di un campione di Bernoulli
# per una griglia di valori del parametro theta (vettore)
# Formula: l(theta) = nr_suc * log(theta) + (n - nr_suc) * log(1 - theta)
# La funzione restituisce l(theta) - max(l(theta)): la curva è traslata
# verso il basso in modo che il massimo cada sempre in 0; questo
# rende confrontabili grafici con numerosità diverse (altrimenti
# la scala verticale cambierebbe completamente al variare di n)
loglik_bern <- function(theta, campione){
  n <- length(campione)
  nr_suc <- sum(campione)
  
  ll <- nr_suc * log(theta) + (n - nr_suc) * log(1 - theta)
  
  return(ll - max(ll))
}


# approssimazione di Taylor del 2° ordine -------------------------------------

# taylor_ll calcola l'approssimazione parabolica della log-verosimiglianza
# ottenuta con lo sviluppo in serie di Taylor arrestato al 2° ordine
# nel punto di massima verosimiglianza theta_hat
#
# Ingredienti:
#   theta_hat  = proporzione campionaria (SMV per Bernoulli)
#   inf_fisher = informazione osservata di Fisher = n / [theta_hat*(1-theta_hat)]
#                (opposto della derivata seconda della log-lik valutata in theta_hat)
#
# Formula: q(theta) = -0.5 * I(theta_hat) * (theta - theta_hat)^2
#
# La parabola ha il vertice in (theta_hat, 0) e si apre verso il basso;
# la curvatura (quanto è "stretta") è determinata dall'informazione di Fisher:
# maggiore è n, più l'informazione è grande e più la parabola è stretta
taylor_ll <- function(theta, campione){
  n <- length(campione)
  # stima di massima verosimiglianza (SMV): media campionaria
  theta_hat <- mean(campione)
  # informazione osservata di Fisher valutata nella SMV
  inf_fisher <- n / (theta_hat * (1 - theta_hat))
  # approssimazione quadratica centrata in theta_hat
  out <- - 0.5 * inf_fisher * (theta - theta_hat)^2
  
  return(out)
}


# grafico esplorativo iniziale ------------------------------------------------

# griglia fitta di valori sullo spazio parametrico [0, 1]
x_theta <- seq(0, 1, by = 0.001)

# campione pilota: n = 10, proporzione campionaria = 0.5
dati <- c(rep(0, 5), rep(1,5))
# log-verosimiglianza normalizzata sulla griglia
y    <- loglik_bern(x_theta, dati)
# approssimazione quadratica di Taylor sulla griglia
tayl <- taylor_ll(x_theta, dati)

# sovrapponiamo i due grafici: curva rossa = log-lik, curva nera = Taylor
plot(x_theta, y,
     type = "l", col = "red", lwd = 3,
     main = "Log-verosimiglianza e approssimazione di Taylor (n = 10, p = 0.5)",
     xlab = expression(theta),
     ylab = expression(l(theta) - l(hat(theta))))
# linea verticale tratteggiata in corrispondenza della SMV (theta_hat)
abline(v = mean(dati), lty = 2)
# sovrappone la parabola di Taylor
lines(x_theta, tayl, lwd = 2)
legend("bottomleft",
       legend = c("log-verosimiglianza", "approx. Taylor 2° ordine", "SMV"),
       col    = c("red", "black", "black"),
       lty    = c(1, 1, 2),
       lwd    = c(3, 2, 1),
       bty    = "n")


# effetto della numerosità n a proporzione fissa (p = 0.5) -------------------
#
# Costruiamo tre campioni con la stessa proporzione campionaria (0.5)
# ma numerosità crescenti (n = 10, 50, 100) per mostrare che
# all'aumentare di n l'informazione di Fisher cresce e la log-lik
# diventa sempre più "a forma di campana": la parabola di Taylor
# diventa un'approssimazione sempre più accurata

dati_10  <- c(rep(0,  5), rep(1,  5))   # n = 10,  theta_hat = 0.5
y_10     <- loglik_bern(x_theta, dati_10)
tayl_10  <- taylor_ll(x_theta,  dati_10)

dati_50  <- c(rep(0, 25), rep(1, 25))   # n = 50,  theta_hat = 0.5
y_50     <- loglik_bern(x_theta, dati_50)
tayl_50  <- taylor_ll(x_theta,  dati_50)

dati_100 <- c(rep(0, 50), rep(1, 50))   # n = 100, theta_hat = 0.5
y_100    <- loglik_bern(x_theta, dati_100)
tayl_100 <- taylor_ll(x_theta,  dati_100)

# griglia 3 × 1: un pannello per ciascuna numerosità
par(mfrow = c(3, 1))

# pannello 1: n = 10
# NOTA: y_10 - max(y_10) è ridondante perché loglik_bern normalizza già
#       internamente (max = 0); il risultato è identico a y_10
plot(x_theta, y_10 - max(y_10),
     type = "l", col = "red", lwd = 3, ylim = c(-20, 0),
     main = "n = 10  |  proporzione campionaria = 0.5",
     xlab = expression(theta),
     ylab = expression(l(theta) - l(hat(theta))))
abline(v = 0.5, lty = 2)
lines(x_theta, tayl_10, lwd = 2)
legend("bottomleft",
       legend = c("log-verosimiglianza", "approx. Taylor 2° ordine", "SMV"),
       col    = c("red", "black", "black"),
       lty    = c(1, 1, 2), lwd = c(3, 2, 1), bty = "n", cex = 0.8)

# pannello 2: n = 50
plot(x_theta, y_50,
     type = "l", col = "red", lwd = 3, ylim = c(-20, 0),
     main = "n = 50  |  proporzione campionaria = 0.5",
     xlab = expression(theta),
     ylab = expression(l(theta) - l(hat(theta))))
abline(v = 0.5, lty = 2)
lines(x_theta, tayl_50, lwd = 2)
legend("bottomleft",
       legend = c("log-verosimiglianza", "approx. Taylor 2° ordine", "SMV"),
       col    = c("red", "black", "black"),
       lty    = c(1, 1, 2), lwd = c(3, 2, 1), bty = "n", cex = 0.8)

# pannello 3: n = 100
plot(x_theta, y_100,
     type = "l", col = "red", lwd = 3, ylim = c(-20, 0),
     main = "n = 100  |  proporzione campionaria = 0.5",
     xlab = expression(theta),
     ylab = expression(l(theta) - l(hat(theta))))
abline(v = 0.5, lty = 2)
lines(x_theta, tayl_100, lwd = 2)
legend("bottomleft",
       legend = c("log-verosimiglianza", "approx. Taylor 2° ordine", "SMV"),
       col    = c("red", "black", "black"),
       lty    = c(1, 1, 2), lwd = c(3, 2, 1), bty = "n", cex = 0.8)

par(mfrow = c(1, 1))


# effetto della numerosità n a proporzione fissa (p = 0.7) -------------------
#
# Ripetiamo l'analisi con proporzione campionaria = 0.7:
# theta_hat si sposta verso destra e la simmetria della parabola
# rimane invariata, ma la curvatura cambia perché l'informazione di Fisher
# dipende da theta_hat*(1-theta_hat)

dati_10  <- c(rep(0,  3), rep(1,  7))   # n = 10,  theta_hat = 0.7
y_10     <- loglik_bern(x_theta, dati_10)
tayl_10  <- taylor_ll(x_theta, dati_10)

dati_50  <- c(rep(0, 15), rep(1, 35))   # n = 50,  theta_hat = 0.7
y_50     <- loglik_bern(x_theta, dati_50)
tayl_50  <- taylor_ll(x_theta, dati_50)

dati_100 <- c(rep(0, 30), rep(1, 70))   # n = 100, theta_hat = 0.7
y_100    <- loglik_bern(x_theta, dati_100)
tayl_100 <- taylor_ll(x_theta, dati_100)

par(mfrow = c(3, 1))

# pannello 1: n = 10
plot(x_theta, y_10 - max(y_10),
     type = "l", col = "red", lwd = 3, ylim = c(-20, 0),
     main = "n = 10  |  proporzione campionaria = 0.7",
     xlab = expression(theta),
     ylab = expression(l(theta) - l(hat(theta))))
abline(v = 0.7, lty = 2)
lines(x_theta, tayl_10, lwd = 2)
legend("bottomleft",
       legend = c("log-verosimiglianza", "approx. Taylor 2° ordine", "SMV"),
       col    = c("red", "black", "black"),
       lty    = c(1, 1, 2), lwd = c(3, 2, 1), bty = "n", cex = 0.8)

# pannello 2: n = 50
plot(x_theta, y_50,
     type = "l", col = "red", lwd = 3, ylim = c(-20, 0),
     main = "n = 50  |  proporzione campionaria = 0.7",
     xlab = expression(theta),
     ylab = expression(l(theta) - l(hat(theta))))
abline(v = 0.7, lty = 2)
lines(x_theta, tayl_50, lwd = 2)
legend("bottomleft",
       legend = c("log-verosimiglianza", "approx. Taylor 2° ordine", "SMV"),
       col    = c("red", "black", "black"),
       lty    = c(1, 1, 2), lwd = c(3, 2, 1), bty = "n", cex = 0.8)

# pannello 3: n = 100
plot(x_theta, y_100,
     type = "l", col = "red", lwd = 3, ylim = c(-20, 0),
     main = "n = 100  |  proporzione campionaria = 0.7",
     xlab = expression(theta),
     ylab = expression(l(theta) - l(hat(theta))))
abline(v = 0.7, lty = 2)
lines(x_theta, tayl_100, lwd = 2)
legend("bottomleft",
       legend = c("log-verosimiglianza", "approx. Taylor 2° ordine", "SMV"),
       col    = c("red", "black", "black"),
       lty    = c(1, 1, 2), lwd = c(3, 2, 1), bty = "n", cex = 0.8)

par(mfrow = c(1, 1))


# effetto della numerosità n a proporzione campionaria = 0.9 -----------------
#
# Quando theta_hat si avvicina ai bordi dello spazio parametrico (0 o 1)
# la log-verosimiglianza diventa sempre più asimmetrica (la coda sinistra
# è compressa verso 0 e quella destra è vincolata al bordo theta = 1);
# la parabola di Taylor, simmetrica per costruzione, fatica quindi
# ad approssimare bene la log-lik con n piccolo
#
dati_10  <- c(rep(0,  1), rep(1,  9))   # n = 10,  theta_hat = 9/10 = 0.9
y_10     <- loglik_bern(x_theta, dati_10)
tayl_10  <- taylor_ll(x_theta, dati_10)

dati_50  <- c(rep(0,  5), rep(1, 45))   # n = 50,  theta_hat = 45/50 = 0.9
y_50     <- loglik_bern(x_theta, dati_50)
tayl_50  <- taylor_ll(x_theta, dati_50)

dati_100 <- c(rep(0, 10), rep(1, 90))   # n = 100, theta_hat = 90/100 = 0.9
y_100    <- loglik_bern(x_theta, dati_100)
tayl_100 <- taylor_ll(x_theta, dati_100)

par(mfrow = c(3, 1))

# pannello 1: n = 10  (theta_hat = 0.9)
# con n piccolo e theta_hat lontano da 0.5 l'asimmetria della log-lik
# è evidente: la parabola non cattura la coda sinistra più allungata
plot(x_theta, y_10 - max(y_10),
     type = "l", col = "red", lwd = 3, ylim = c(-20, 0),
     main = "n = 10  |  proporzione campionaria = 0.9",
     xlab = expression(theta),
     ylab = expression(l(theta) - l(hat(theta))))
abline(v = 0.9, lty = 2)
lines(x_theta, tayl_10, lwd = 2)
legend("bottomleft",
       legend = c("log-verosimiglianza", "approx. Taylor 2° ordine", "SMV"),
       col    = c("red", "black", "black"),
       lty    = c(1, 1, 2), lwd = c(3, 2, 1), bty = "n", cex = 0.8)

# pannello 2: n = 50
# con n più grande la log-lik si assottiglia e diventa più simmetrica:
# l'approssimazione parabolica migliora visibilmente
plot(x_theta, y_50,
     type = "l", col = "red", lwd = 3, ylim = c(-20, 0),
     main = "n = 50  |  proporzione campionaria = 0.9",
     xlab = expression(theta),
     ylab = expression(l(theta) - l(hat(theta))))
abline(v = 0.9, lty = 2)
lines(x_theta, tayl_50, lwd = 2)
legend("bottomleft",
       legend = c("log-verosimiglianza", "approx. Taylor 2° ordine", "SMV"),
       col    = c("red", "black", "black"),
       lty    = c(1, 1, 2), lwd = c(3, 2, 1), bty = "n", cex = 0.8)

# pannello 3: n = 100
# con n grande l'approssimazione quadratica è praticamente sovrapposta
# alla log-verosimiglianza: il teorema asintotico garantisce che la
# log-lik standardizzata converge ad una forma normale / parabolica
plot(x_theta, y_100,
     type = "l", col = "red", lwd = 3, ylim = c(-20, 0),
     main = "n = 100  |  proporzione campionaria = 0.9",
     xlab = expression(theta),
     ylab = expression(l(theta) - l(hat(theta))))
abline(v = 0.9, lty = 2)
lines(x_theta, tayl_100, lwd = 2)
legend("bottomleft",
       legend = c("log-verosimiglianza", "approx. Taylor 2° ordine", "SMV"),
       col    = c("red", "black", "black"),
       lty    = c(1, 1, 2), lwd = c(3, 2, 1), bty = "n", cex = 0.8)

par(mfrow = c(1, 1))