# iterazione con map2, pmap e pwalk (pacchetto purrr) ------------------------
library(purrr)


# map2: iterare su due vettori in parallelo -----------------------------------

# map2 applica una funzione a due vettori elemento per elemento:
# il primo elemento di .x è abbinato al primo di .y, il secondo al secondo, ...
# .f = ~c(.x, .y) crea semplicemente la coppia di valori: utile per
# verificare visivamente come vengono abbinati gli argomenti
medie <- c(10, 20, 30)
sqm   <- c(1, 2, 3)

map2(medie, sqm, .f = ~c(.x, .y))

# applicazione pratica: genero 3 campioni normali con media e deviazione
# standard diverse; l'output è una lista di 3 vettori da 100 elementi ciascuno
out <- map2(.x = medie, .y = sqm, .f = ~rnorm(n = 100, mean = .x, sd = .y))

# map_vec semplifica l'output da lista a vettore atomico:
# calcolo la media e la deviazione standard di ciascun campione generato
map_vec(out, mean)
map_vec(out, sd)


# pmap: iterare su tre o più vettori in parallelo -----------------------------

# quando gli argomenti da variare sono più di due si usa pmap:
# .l accetta una lista (o un data.frame) di vettori della stessa lunghezza;
# dentro .f gli argomenti si richiamano con ..1, ..2, ..3 in base
# alla posizione nella lista (o con il nome della colonna del data.frame)
n <- c(5, 50, 100)

# input come lista: ogni elemento della lista è un vettore di argomenti
list(n = n, medie = medie, sqm = sqm)
pmap(.l = list(n = n, medie = medie, sqm = sqm), 
      .f = ~rnorm(..1, ..2, ..3))

# input come data.frame: equivalente alla lista ma più leggibile
# quando gli argomenti hanno una struttura tabulare naturale
data.frame(n = n, medie = medie, sqm = sqm)
pmap(.l = data.frame(n = n, medie = medie, sqm = sqm), 
     .f = ~rnorm(..1, ..2, ..3))


# pwalk: come pmap ma per effetti collaterali (grafici, stampe, salvataggi) ---

# walk (e le sue varianti walk2, pwalk) è la versione di map pensata per
# le funzioni che vengono chiamate per il loro effetto collaterale e non
# per il valore restituito; in questo caso l'effetto è disegnare un istogramma
# NOTA: pwalk non restituisce nulla di utile (.l viene restituito
#       invisibilmente): lo usiamo solo per i grafici che produce
list(n = n, medie = medie, sqm = sqm)
pwalk(.l = data.frame(n = n, medie = medie, sqm = sqm), 
      .f = ~hist(rnorm(..1, ..2, ..3)))


# funzioni dalla lezione precedente -------------------------------------------

# log-verosimiglianza normalizzata di Bernoulli su una griglia di valori theta
# (sottraendo il massimo la curva picca sempre in 0: scale confrontabili)
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 nel punto di massima verosimiglianza:
# parabola con vertice in theta_hat = 0 e curvatura pari all'informazione
# osservata di Fisher I(theta_hat) = n / [theta_hat * (1 - theta_hat)]
taylor_ll <- function(theta, campione){
  n         <- length(campione)
  theta_hat <- mean(campione)
  inf_fisher <- n / (theta_hat * (1 - theta_hat))
  
  out <- -0.5 * inf_fisher * (theta - theta_hat)^2
  
  return(out)
}

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

# calcolo su un campione pilota (n = 10, p = 0.5) per verificare le funzioni
dati <- c(rep(0, 5), rep(1, 5))
y    <- loglik_bern(x_theta, dati)
tayl <- taylor_ll(x_theta, dati)

# grafico di controllo: le due curve dovrebbero sovrapporsi bene vicino al picco
plot(x_theta, y, type = "l", col = "red", lwd = 3)
abline(v = mean(dati))
lines(x_theta, tayl)


# genera_campione: funzione ausiliaria ----------------------------------------

# genera_campione costruisce un vettore di 0 e 1 con la proporzione di
# successi (pi_hat) desiderata e la numerosità n indicata;
# round() garantisce che il numero di successi sia un intero
# (es. n = 10, pi_hat = 0.7  -->  7 uni e 3 zeri)
genera_campione <- function(n, pi_hat){
  n_suc <- round(n * pi_hat)
  dati  <- c(rep(0, n - n_suc), rep(1, n_suc))
  
  return(dati)
}


# riscrittura con programmazione funzionale -----------------------------------
#
# Nella lezione precedente i tre pannelli erano ottenuti con tre blocchi
# di codice identici (uno per n = 10, uno per n = 50, uno per n = 100).
# Con map e walk2 si elimina la ripetizione: si definisce la logica
# una sola volta e la si applica a tutti i valori di sample_size.

sample_size <- c(10, 50, 100)

# passo 1: genero i tre campioni con map (output: lista di 3 vettori)
dati <- map(sample_size, genera_campione, pi_hat = 0.5)
# passo 2: calcolo log-verosimiglianza e Taylor per ciascun campione
ll_bernoulli   <- map(dati, ~loglik_bern(theta = x_theta, campione = .))
tayl_bernoulli <- map(dati, ~taylor_ll(theta = x_theta,   campione = .))

# NOTA: ylim = c(-2, 0) è uno zoom attorno al picco della curva;
#       con numerosità piccole la log-lik scende molto più in basso,
#       ma qui ci interessa vedere il comportamento vicino alla SMV
# NOTA: lo stesso ylim per tutti i pannelli rende il confronto tra
#       le numerosità diretto, a differenza di lasciare la scala libera

# --- percorso incrementale: come si costruisce un grafico con R base --------


# passo 1 - grafico minimo (nessun titolo, nessuna legenda) -------------------
# plot() con type = "l" disegna la curva; lines() sovrappone una seconda curva
# allo stesso grafico senza aprirne uno nuovo
par(mfrow = c(3, 1))
walk2(.x = ll_bernoulli,
      .y = tayl_bernoulli,
      .f = ~{
        plot(x_theta, .x, type = "l", col = "red", lwd = 3, ylim = c(-2, 0))
        lines(x_theta, .y)
      })
par(mfrow = c(1, 1))


# passo 2 - aggiunta di titolo e etichette degli assi ------------------------
# main: titolo del pannello
# xlab, ylab: etichette degli assi; expression() permette di usare
# simboli matematici (theta, l(theta) con il cappello)
par(mfrow = c(3, 1))
walk2(.x = ll_bernoulli,
      .y = tayl_bernoulli,
      .f = ~{
        plot(x_theta, .x,
             type = "l", col = "red", lwd = 3, ylim = c(-2, 0),
             main = "Log-verosimiglianza e approssimazione di Taylor",
             xlab = expression(theta),
             ylab = expression(l(theta) - l(hat(theta))))
        lines(x_theta, .y)
      })
par(mfrow = c(1, 1))


# passo 3 - linea verticale per la SMV con stile personalizzato ---------------
# abline(v = ...) aggiunge una retta verticale; lty = 2 la rende tratteggiata
# per distinguerla visivamente dalle curve
# mean(.x) non funziona qui perché .x è la log-lik già calcolata (non il campione);
# il massimo della log-lik normalizzata è 0: x_theta[which.max(.x)] individua
# il valore di theta a cui corrisponde il massimo
par(mfrow = c(3, 1))
walk2(.x = ll_bernoulli,
      .y = tayl_bernoulli,
      .f = ~{
        plot(x_theta, .x,
             type = "l", col = "red", lwd = 3, ylim = c(-2, 0),
             main = "Log-verosimiglianza e approssimazione di Taylor",
             xlab = expression(theta),
             ylab = expression(l(theta) - l(hat(theta))))
        lines(x_theta, .y)
        abline(v = x_theta[which.max(.x)], lty = 2)
      })
par(mfrow = c(1, 1))


# passo 4 - aggiunta della legenda --------------------------------------------
# legend() posiziona la legenda in un angolo del grafico:
#   - primo argomento: posizione ("bottomleft", "topright", ecc.)
#   - legend: etichette testuali
#   - col, lty, lwd: devono corrispondere agli stessi argomenti usati
#     in plot(), lines() e abline()
#   - bty = "n": elimina il bordo rettangolare attorno alla legenda
#   - cex = 0.8: riduce la dimensione del testo per non occupare troppo spazio
par(mfrow = c(3, 1))
walk2(.x = ll_bernoulli,
      .y = tayl_bernoulli,
      .f = ~{
        plot(x_theta, .x,
             type = "l", col = "red", lwd = 3, ylim = c(-2, 0),
             main = "Log-verosimiglianza e approssimazione di Taylor",
             xlab = expression(theta),
             ylab = expression(l(theta) - l(hat(theta))))
        lines(x_theta, .y, lwd = 2)
        abline(v = x_theta[which.max(.x)], lty = 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))


# passo 5 - titolo personalizzato per pannello con pwalk ----------------------
# Per includere n nel titolo di ciascun pannello dobbiamo passare un terzo
# argomento variabile; walk2 ne gestisce solo due (.x e .y), quindi
# passiamo a pwalk che accetta una lista di lunghezza arbitraria.
# ..3 (o equivalentemente il nome "n" assegnato nella lista) contiene
# la numerosità corrispondente a ciascuna iterazione
par(mfrow = c(3, 1))
pwalk(.l = list(ll   = ll_bernoulli,
                tayl = tayl_bernoulli,
                n    = sample_size),
      .f = ~{
        plot(x_theta, ..1,
             type = "l", col = "red", lwd = 3, ylim = c(-2, 0),
             main = paste0("n = ", ..3, "  |  proporzione campionaria = 0.5"),
             xlab = expression(theta),
             ylab = expression(l(theta) - l(hat(theta))))
        lines(x_theta, ..2, lwd = 2)
        abline(v = x_theta[which.max(..1)], lty = 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))


# stesso schema con proporzione campionaria = 0.9 ----------------------------
# cambia solo pi_hat in genera_campione: il resto del codice è identico,
# e questo è il vantaggio della programmazione funzionale rispetto alla
# ripetizione manuale dei blocchi

dati           <- map(sample_size, genera_campione, pi_hat = 0.9)
ll_bernoulli   <- map(dati, ~loglik_bern(theta = x_theta, campione = .))
tayl_bernoulli <- map(dati, ~taylor_ll(theta = x_theta,   campione = .))

par(mfrow = c(3, 1))
pwalk(.l = list(ll   = ll_bernoulli,
                tayl = tayl_bernoulli,
                n    = sample_size),
      .f = ~{
        plot(x_theta, ..1,
             type = "l", col = "red", lwd = 3, ylim = c(-2, 0),
             main = paste0("n = ", ..3, "  |  proporzione campionaria = 0.9"),
             xlab = expression(theta),
             ylab = expression(l(theta) - l(hat(theta))))
        lines(x_theta, ..2, lwd = 2)
        abline(v = x_theta[which.max(..1)], lty = 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))