# verosimiglianza, log-verosimiglianza e verosimiglianza relativa
# per il modello Uniforme(0, theta) ----------------------------------------
#
# Il parametro theta è il limite superiore dell'intervallo: la SMV è
# il massimo campionario x_(n), perché è il più piccolo valore di theta
# compatibile con i dati osservati.


# funzione di verosimiglianza -------------------------------------------------

# L(theta) = theta^(-n)  se theta >= x_(n)   (tutti i dati "stanno dentro")
# L(theta) = 0           se theta < x_(n)    (almeno un dato supera theta)
# la funzione è definita su una griglia di valori theta (vettore):
# ifelse gestisce la condizione elemento per elemento
lik_unif <- function(theta, campione) {
  n     <- length(campione)
  x_max <- max(campione)
  
  l <- ifelse(theta < x_max, 0, theta^(-n))
  return(l)
}

# campione simulato di n = 10 osservazioni da Uniforme(0, 5)
# set.seed garantisce la riproducibilità del campione casuale
set.seed(42)
dati  <- runif(10, min = 0, max = 5)
# massimo campionario: è la SMV di theta
x_max <- max(dati)

# griglia di theta su cui valutare la verosimiglianza
# si parte da 0 per mostrare il salto in corrispondenza di x_max
x_theta <- seq(0, 10, by = 0.001)
y       <- lik_unif(x_theta, dati)

# grafico della verosimiglianza: si vede chiaramente il salto da 0 a
# theta^(-n) in corrispondenza di x_max, e la successiva discesa monotona
plot(x_theta, y,
     type = "l", col = "red", lwd = 3,
     xlab = expression(theta),
     ylab = expression(L(theta)),
     main = "Verosimiglianza - Uniforme(0, theta)")

# linea verticale tratteggiata in corrispondenza della SMV (massimo campionario)
abline(v = x_max, col = "blue",      lty = 2, lwd = 2)
# linea verticale in corrispondenza del valore vero del parametro
abline(v = 5,     col = "darkgreen", lty = 2, lwd = 2)

# expression() nella legenda permette di usare notazione matematica:
# hat(theta) produce theta con il cappello; x[(n)] produce x con pedice (n)
legend("topright",
       legend = c(expression(hat(theta) == x[(n)]),
                  expression(theta[vero] == 5)),
       col = c("blue", "darkgreen"),
       lty = 2, lwd = 2)


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

# l(theta) = -n * log(theta)  se theta >= x_(n)
# l(theta) = -Inf              se theta < x_(n)   (log(0) = -Inf)
# -Inf è il valore corretto dal punto di vista matematico e consente
# di usare ifelse senza valori speciali
loglik_unif <- function(theta, campione) {
  n     <- length(campione)
  x_max <- max(campione)
  
  ll <- ifelse(theta < x_max, -Inf, -n * log(theta))
  return(ll)
}

y_log <- loglik_unif(x_theta, dati)

# ylim = c(-20, 0) taglia i valori -Inf a sinistra di x_max:
# senza questo limite il grafico sarebbe illeggibile
plot(x_theta, y_log,
     type = "l", col = "red", lwd = 3,
     ylim = c(-20, 0),
     xlab = expression(theta),
     ylab = expression(l(theta)),
     main = "Log-verosimiglianza - Uniforme(0, theta)")

abline(v = x_max, col = "blue",      lty = 2, lwd = 2)
abline(v = 5,     col = "darkgreen", lty = 2, lwd = 2)

legend("topright",
       legend = c(expression(hat(theta) == x[(n)]),
                  expression(theta[vero] == 5)),
       col = c("blue", "darkgreen"),
       lty = 2, lwd = 2)


# verosimiglianza relativa e log-verosimiglianza relativa ---------------------

# R(theta) = L(theta) / L(theta_hat) = (x_(n) / theta)^n  se theta >= x_(n)
# è normalizzata in [0, 1] e vale 1 nel punto di massima verosimiglianza;
# è utile per confrontare campioni con numerosità diverse perché
# la scala assoluta di L(theta) dipende da n
rel_lik_unif <- function(theta, campione) {
  n     <- length(campione)
  x_max <- max(campione)
  
  R <- ifelse(theta < x_max, 0, (x_max / theta)^n)
  return(R)
}

# r(theta) = log R(theta) = n * log(x_(n) / theta)  se theta >= x_(n)
# equivalente alla log-verosimiglianza traslata in modo che il massimo sia 0;
# stesso risultato che si ottiene con loglik_bern(...) - max(loglik_bern(...))
# nel caso della Bernoulli
rel_loglik_unif <- function(theta, campione) {
  n     <- length(campione)
  x_max <- max(campione)
  
  r <- ifelse(theta < x_max, -Inf, n * log(x_max / theta))
  return(r)
}


# confronto tra numerosità diverse sullo stesso grafico ----------------------
#
# Sovrapporre più curve sullo stesso pannello richiede di:
# 1) aprire il grafico vuoto con plot(NULL, ...) impostando gli assi
# 2) aggiungere una curva per volta con lines()
# Questo schema si presta naturalmente all'iterazione con un ciclo for
# o con le funzioni walk della famiglia purrr

theta_vero <- 5
numerosita <- c(5, 10, 30, 100)
set.seed(42)

x_theta <- seq(0, 10, by = 0.001)
# un colore distinto per ciascuna numerosità: stesso ordine di numerosita
colori  <- c("red", "blue", "darkgreen", "purple")


# --- verosimiglianza relativa: ciclo for ------------------------------------

# plot(NULL, ...) apre un piano cartesiano vuoto con assi e titolo pronti;
# le curve vengono aggiunte dall'iterazione successiva
plot(NULL,
     xlim = c(4, 10), ylim = c(0, 1),
     xlab = expression(theta),
     ylab = expression(R(theta)),
     main = expression("Verosimiglianza relativa - Uniforme(0," ~ theta ~ ")"))

# seq_along(numerosita) genera l'indice k = 1, 2, 3, 4;
# usare l'indice permette di accedere contemporaneamente a numerosita[k]
# e a colori[k] con la stessa posizione
# NOTA: i campioni non sono riproducibili fra un'esecuzione e l'altra
#       perché set.seed(42) è stato chiamato una sola volta prima del ciclo;
#       ogni iterazione consuma il generatore casuale in sequenza
for (k in seq_along(numerosita)) {
  n    <- numerosita[k]
  dati <- runif(n, 0, theta_vero)
  y    <- rel_lik_unif(x_theta, dati)
  lines(x_theta, y, col = colori[k], lwd = 2)
}

# linea verticale al valore vero del parametro (riferimento fisso)
abline(v = theta_vero, lty = 2, col = "black", lwd = 1.5)
# paste0("n = ", numerosita) costruisce automaticamente le etichette
# "n = 5", "n = 10", "n = 30", "n = 100"
legend("topright",
       legend = paste0("n = ", numerosita),
       col = colori, lwd = 2)


# --- verosimiglianza relativa: variante con walk2 (purrr) -------------------
#
# walk2 itera su due vettori in parallelo (numerosita e colori) ed è
# l'equivalente funzionale del ciclo for con seq_along;
# .x riceve la numerosità, .y il colore corrispondente
# NOTA: set.seed(42) deve essere impostato di nuovo prima di walk2 per
#       ottenere gli stessi campioni del ciclo for

plot(NULL,
     xlim = c(4, 10), ylim = c(0, 1),
     xlab = expression(theta),
     ylab = expression(R(theta)),
     main = expression("Verosimiglianza relativa - Uniforme(0," ~ theta ~ ")"))

library(purrr)
set.seed(42)
walk2(numerosita, colori,
      ~{
        dati <- runif(.x, 0, theta_vero)
        y    <- rel_lik_unif(x_theta, dati)
        lines(x_theta, y, col = .y, lwd = 2)
      })

abline(v = theta_vero, lty = 2, col = "black", lwd = 1.5)
legend("topright",
       legend = paste0("n = ", numerosita),
       col = colori, lwd = 2)


# --- log-verosimiglianza relativa: ciclo for ---------------------------------

plot(NULL,
     xlim = c(4, 10), ylim = c(-10, 0),
     xlab = expression(theta),
     ylab = expression(r(theta)),
     main = expression("Log-verosimiglianza relativa - Uniforme(0," ~ theta ~ ")"))

# in questo ciclo si usa set.seed(k) ad ogni iterazione: ogni numerosità
# riceve un campione casuale indipendente e riproducibile (seed 1, 2, 3, 4)
for (k in seq_along(numerosita)) {
  n <- numerosita[k]
  set.seed(k)
  dati <- runif(n, 0, theta_vero)
  y    <- rel_loglik_unif(x_theta, dati)
  lines(x_theta, y, col = colori[k], lwd = 2)
}

abline(v = theta_vero, lty = 2, col = "black", lwd = 1.5)
# abline(h = 0) aggiunge una retta orizzontale in corrispondenza di r = 0,
# ovvero il livello del massimo; aiuta a leggere la curvatura delle curve
abline(h = 0, lty = 3, col = "gray50")
legend("topright",
       legend = paste0("n = ", numerosita),
       col = colori, lwd = 2)


# --- log-verosimiglianza relativa: variante con pwalk (purrr) ---------------
#
# Qui servono tre argomenti variabili (n, colore, seed): walk2 non basta,
# si usa pwalk con una lista a tre elementi; dentro .f gli argomenti
# si richiamano con ..1 (n), ..2 (colore), ..3 (seed)

plot(NULL,
     xlim = c(4, 10), ylim = c(-10, 0),
     xlab = expression(theta),
     ylab = expression(r(theta)),
     main = expression("Log-verosimiglianza relativa - Uniforme(0," ~ theta ~ ")"))

pwalk(list(n     = numerosita,
           col   = colori,
           seed  = seq_along(numerosita)),
      ~{
        set.seed(..3)
        dati <- runif(..1, 0, theta_vero)
        y    <- rel_loglik_unif(x_theta, dati)
        lines(x_theta, y, col = ..2, lwd = 2)
      })

abline(v = theta_vero, lty = 2, col = "black", lwd = 1.5)
abline(h = 0, lty = 3, col = "gray50")
legend("topright",
       legend = paste0("n = ", numerosita),
       col = colori, lwd = 2)


# campioni con massimo campionario fissato ------------------------------------
#
# Fissare x_(n) permette di isolare l'effetto di n sulla curvatura della
# verosimiglianza relativa, eliminando la variabilità dovuta al fatto che
# campioni più grandi tendono ad avere massimi campionari più vicini a theta_vero

theta_vero  <- 5
x_max_fisso <- 4.8   # massimo campionario identico per tutti i campioni
numerosita  <- c(5, 10, 30, 100)
set.seed(42)

# genera_campione costruisce un campione di n osservazioni con massimo esatto
# pari a x_max: genera n-1 valori casuali in (0, x_max) e aggiunge x_max
# NOTA: il parametro x_max della funzione non viene usato internamente;
#       la funzione fa riferimento alla variabile globale x_max_fisso.
#       Nella pratica è preferibile usare il parametro formale per evitare
#       dipendenze da variabili esterne (effetti collaterali nascosti)
genera_campione <- function(n, x_max, theta_vero) {
  c(runif(n - 1, min = 0, max = x_max_fisso), x_max)
}

x_theta <- seq(0, 10, by = 0.001)
colori  <- c("red", "blue", "darkgreen", "purple")


# --- verosimiglianza relativa a massimo fissato: ciclo for ------------------

plot(NULL,
     xlim = c(4, 10), ylim = c(0, 1),
     xlab = expression(theta),
     ylab = expression(R(theta)),
     main = expression("Verosimiglianza relativa - Uniforme(0," ~ theta ~ ")"))

for (k in seq_along(numerosita)) {
  dati <- genera_campione(numerosita[k], x_max_fisso, theta_vero)
  y    <- rel_lik_unif(x_theta, dati)
  lines(x_theta, y, col = colori[k], lwd = 2)
}

abline(v = theta_vero,  lty = 2, col = "black",  lwd = 1.5)
abline(v = x_max_fisso, lty = 3, col = "gray40", lwd = 1.5)
# la legenda include sia le curve per numerosità sia le due linee verticali
# (theta_vero e x_(n)): le due voci extra richiedono di specificare
# lty separatamente per ciascun elemento (1 per le curve, 2 e 3 per le linee)
legend("topright",
       legend = c(paste0("n = ", numerosita),
                  expression(theta[vero]), expression(x[(n)])),
       col    = c(colori, "black", "gray40"),
       lty    = c(1, 1, 1, 1, 2, 3),
       lwd    = 2)


# --- verosimiglianza relativa a massimo fissato: variante con walk2 ---------

plot(NULL,
     xlim = c(4, 10), ylim = c(0, 1),
     xlab = expression(theta),
     ylab = expression(R(theta)),
     main = expression("Verosimiglianza relativa - Uniforme(0," ~ theta ~ ")"))

walk2(numerosita, colori,
      ~{
        dati <- genera_campione(.x, x_max_fisso, theta_vero)
        y    <- rel_lik_unif(x_theta, dati)
        lines(x_theta, y, col = .y, lwd = 2)
      })

abline(v = theta_vero,  lty = 2, col = "black",  lwd = 1.5)
abline(v = x_max_fisso, lty = 3, col = "gray40", lwd = 1.5)
legend("topright",
       legend = c(paste0("n = ", numerosita),
                  expression(theta[vero]), expression(x[(n)])),
       col    = c(colori, "black", "gray40"),
       lty    = c(1, 1, 1, 1, 2, 3),
       lwd    = 2)


# --- log-verosimiglianza relativa a massimo fissato: ciclo for --------------

plot(NULL,
     xlim = c(4, 10), ylim = c(-10, 0),
     xlab = expression(theta),
     ylab = expression(r(theta)),
     main = expression("Log-verosimiglianza relativa - Uniforme(0," ~ theta ~ ")"))

for (k in seq_along(numerosita)) {
  dati <- genera_campione(numerosita[k], x_max_fisso, theta_vero)
  y    <- rel_loglik_unif(x_theta, dati)
  lines(x_theta, y, col = colori[k], lwd = 2)
}

abline(v = theta_vero,  lty = 2, col = "black",  lwd = 1.5)
# nota: il colore della linea per x_(n) è "orange" nel codice originale
# ma "gray40" nella legenda; qui si mantiene la coerenza usando "gray40"
abline(v = x_max_fisso, lty = 3, col = "gray40", lwd = 1.5)
abline(h = 0,           lty = 3, col = "gray50")
legend("topright",
       legend = c(paste0("n = ", numerosita),
                  expression(theta[vero]), expression(x[(n)])),
       col    = c(colori, "black", "gray40"),
       lty    = c(1, 1, 1, 1, 2, 3),
       lwd    = 2)


# --- log-verosimiglianza relativa a massimo fissato: variante con walk2 -----

plot(NULL,
     xlim = c(4, 10), ylim = c(-10, 0),
     xlab = expression(theta),
     ylab = expression(r(theta)),
     main = expression("Log-verosimiglianza relativa - Uniforme(0," ~ theta ~ ")"))

walk2(numerosita, colori,
      ~{
        dati <- genera_campione(.x, x_max_fisso, theta_vero)
        y    <- rel_loglik_unif(x_theta, dati)
        lines(x_theta, y, col = .y, lwd = 2)
      })

abline(v = theta_vero,  lty = 2, col = "black",  lwd = 1.5)
abline(v = x_max_fisso, lty = 3, col = "gray40", lwd = 1.5)
abline(h = 0,           lty = 3, col = "gray50")
legend("topright",
       legend = c(paste0("n = ", numerosita),
                  expression(theta[vero]), expression(x[(n)])),
       col    = c(colori, "black", "gray40"),
       lty    = c(1, 1, 1, 1, 2, 3),
       lwd    = 2)