leapfrog_ext <- function(Z, beta, mu, gamma, P_Z, P_beta, P_mu, P_gamma, Y, X, site_idx, eps, L, f_grad_Z, f_grad_beta, f_grad_mu, f_grad_gamma) {
# Premier demi-pas pour les impulsions (Momentum half-step)
P_Z <- P_Z - 0.5 * eps * f_grad_Z(Z, beta, Y, X, site_idx, mu, gamma)
P_beta <- P_beta - 0.5 * eps * f_grad_beta(Z, beta, Y, X, site_idx)
P_mu <- P_mu - 0.5 * eps * f_grad_mu(Z, mu, gamma)
P_gamma <- P_gamma - 0.5 * eps * f_grad_gamma(Z, mu, gamma)
for (l in 1:L) {
# Mise à jour des positions (Position update)
Z <- Z + eps * P_Z
beta <- beta + eps * P_beta
mu <- mu + eps * P_mu
gamma <- gamma + eps * P_gamma
if (l < L) {
# Mise à jour complète des impulsions (sauf pour le dernier pas)
P_Z <- P_Z - eps * f_grad_Z(Z, beta, Y, X, site_idx, mu, gamma)
P_beta <- P_beta - eps * f_grad_beta(Z, beta, Y, X, site_idx)
P_mu <- P_mu - eps * f_grad_mu(Z, mu, gamma)
P_gamma <- P_gamma - eps * f_grad_gamma(Z, mu, gamma)
}
}
# Dernier demi-pas pour les impulsions (Final momentum half-step)
P_Z <- P_Z - 0.5 * eps * f_grad_Z(Z, beta, Y, X, site_idx, mu, gamma)
P_beta <- P_beta - 0.5 * eps * f_grad_beta(Z, beta, Y, X, site_idx)
P_mu <- P_mu - 0.5 * eps * f_grad_mu(Z, mu, gamma)
P_gamma <- P_gamma - 0.5 * eps * f_grad_gamma(Z, mu, gamma)
return(list(Z = Z, beta = beta, mu = mu, gamma = gamma, P_Z = P_Z, P_beta = P_beta, P_mu = P_mu, P_gamma = P_gamma))
}
HMC_step_MC_ext <- function(Z, beta, mu, gamma, Y, X, site_idx, eps, L, lp_func, f_grad_Z, f_grad_beta, f_grad_mu, f_grad_gamma) {
# 1. Initialisation des impulsions aleatoire
P_Z_init <- rnorm(length(Z))
P_beta_init <- rnorm(1)
P_mu_init <- rnorm(1)
P_gamma_init <- rnorm(1)
# 2. Calcul de l'énergie initiale (Hamiltonien initial)
init_U <- -lp_func(Z, beta, Y, X, site_idx, mu, gamma)
init_K <- sum(P_Z_init^2)/2 + P_beta_init^2/2 + P_mu_init^2/2 + P_gamma_init^2/2
# 3. Leapfrog
res <- leapfrog_ext(Z, beta, mu, gamma, P_Z_init, P_beta_init, P_mu_init, P_gamma_init, Y, X, site_idx, eps, L, f_grad_Z, f_grad_beta, f_grad_mu, f_grad_gamma)
# 4. Calcul de l'énergie finale (Hamiltonien final)
final_U <- -lp_func(res$Z, res$beta, Y, X, site_idx, res$mu, res$gamma)
final_K <- sum(res$P_Z^2)/2 + res$P_beta^2/2 + res$P_mu^2/2 + res$P_gamma^2/2
# 5. Accept/Reject
if (log(runif(1)) < (init_U + init_K - final_U - final_K)) {
return(list(Z = res$Z, beta = res$beta, mu = res$mu, gamma = res$gamma, accepted = TRUE))
} else {
return(list(Z = Z, beta = beta, mu = mu, gamma = gamma, accepted = FALSE))
}
}
run_chain_ext <- function(id, n_iter, Y, X, site_idx, eps, L) {
# 1. Initialisation Empirique (Pour favoriser la convergence)
# Calcul de la moyenne par site pour mồi (moyenner sur les années)
# Correction : on utilise Y (l'argument de la fonction) et non Y_ext (variable globale)
Y_mean_site <- as.numeric(tapply(Y, site_idx, mean))
# Initialisation de Z (log-moyenne par site)
curr_Z <- log(Y_mean_site + 0.1)
# Initialisation de beta (effet neutre au départ)
curr_beta <- 0
# Initialisation des hyperparamètres (mu et gamma)
curr_mu <- mean(curr_Z)
curr_gamma <- log(var(curr_Z))
# 2. Sauvegarde des échantillons (Traces)
beta_trace <- numeric(n_iter)
mu_trace <- numeric(n_iter)
sigma_trace <- numeric(n_iter)
Z_trace <- matrix(0, n_iter, 23)
accepted_count_ext <- 0
# 3. Boucle MCMC (Algorithme HMC)
for (k in 1:n_iter) {
res <- HMC_step_MC_ext(curr_Z, curr_beta, curr_mu, curr_gamma, Y, X, site_idx, eps, L,
lp_func = log_posterior_ext,
f_grad_Z = grad_Z_ext,
f_grad_beta = grad_beta,
f_grad_mu = grad_mu_ext,
f_grad_gamma = grad_gamma_ext)
# Mise à jour des valeurs actuelles
curr_Z <- res$Z
curr_beta <- res$beta
curr_mu <- res$mu
curr_gamma <- res$gamma
if (res$accepted) {
accepted_count_ext <- accepted_count_ext + 1
}
# Enregistrement dans les traces
beta_trace[k] <- curr_beta
mu_trace[k] <- curr_mu
sigma_trace[k] <- sqrt(exp(curr_gamma))
Z_trace[k, ] <- curr_Z
# Affichage de la progression et du taux d'acceptation
if (k %% max(1, floor(n_iter / 10)) == 0) {
current_acc_rate <- (accepted_count_ext / k) * 100
cat(sprintf(" Chaîne %d | Progression : %d / %d (%d%%) | Taux d'acc. : %.1f%%\n",
id, k, n_iter, round(k/n_iter*100), current_acc_rate))
}
}
final_acc_rate_ext <- (accepted_count_ext / n_iter)
cat(sprintf("=> Chaîne %d terminée avec succès. Taux d'acceptation final : %.2f%%\n",
id, final_acc_rate_ext * 100))
return(list(mu = mu_trace, sigma = sigma_trace, beta = beta_trace, Z = Z_trace, acc_rate = final_acc_rate_ext))
}
