---
title: "MCMC"
format: html
editor: visual
---

## Premier exemple de Metropolis-Hastings

On va générer un échantillon selon une $Exp(\lambda)$ avec $\lambda = 3$ en utilisant l'algorithme de Metropolis–Hastings (M.-H.)

```{r}
nchain <- 10000## longueur de la chaine de Markov
chain <- rep(NA, nchain)## la chaine
init <- 10## etat initial de la chaine
sigma <- 1## ecart-type du noyau de proposition (Gaussien)
lambda <- 3## paramètre de la loi Exponentielle que l'on souhaite simuler

chain[1] <- init
for (i in 1:(nchain-1)){
  prop <- rnorm(1, chain[i], sigma)##candidat
  
  tmp <- exp(dexp(prop, lambda, log = TRUE) -
               dexp(chain[i], lambda, log = TRUE))##je calcule via les log-densités car plus robuste numériquement
  acc_prob <- min(1, tmp)## faites moi confiance j'ai simplifié les maths de l'algo
  
  if (runif(1) < acc_prob){
    chain[i+1] <- prop
  } else
    chain[i+1] <- chain[i]
}

plot(chain, type = "l", xlab = "Iteration",
     ylab = expression(X[t]))
```

On voit que la chaine a atteint son régime stationnaire après 3000 itérations, regardons cela.

```{r}
library(coda)##à installer d'abord
plot(mcmc(chain, start = 3000))##je vire les 3000 1ère observations de la chaine
```

On va comparer l'histogramme de ma chaine de Markov avec la densité théorique $Exp(\lambda)$.

```{r}
hist(chain[-(1:3000)], freq = FALSE)##je vire les 3000 1ère obs
g <- function(x) dexp(x, lambda)
curve(g, from = 0, to = 12, col = 2, add = TRUE)

```

## M.-H. avec les extrêmes :

Ce qu'on va faire c'est ajuster une GEV dans un cadre Bayésien. Plus précisement, on va s'intéresser au modèle suivant

$$
M_i \mid \{\mu, \sigma, \xi\} \stackrel{iid}{\sim} GEV(\mu, \sigma, \xi),
$$

et l'on prendra comme loi a priori pour \$\\mu, \sigma\$ et $\xi$ des loi normales indépendantes (sauf pour $\sigma$ où se sera une log-normale car il est positif).

Pour cela on va utiliser la library **evdbayes** (que l'on installe d'abord)

```{r}
library(evdbayes)
Sigma <- diag(c(10000, 10000, 100))
prior <- prior.norm(mean = c(0,0,0), cov = Sigma)
```

Les lignes de code plus haut suppose que l'on a a priori des loi normales (ou log-normal pour $\sigma$) indépendantes, centrées et de variance 10000, 10000 et 100.

```{r}
url <- "https://mribatet.perso.math.cnrs.fr/docs/EVT/TG_SOUID260218.txt"
data <- read.table(url, skip = 20,
                   header = TRUE,
                   sep = ",", 
                   na.strings = "-9999")
data$TG <- data$TG / 10##pour être en degrés Celsius

data$DATE <- as.Date(as.character(data$DATE), 
                     format = "%Y%m%d") 

library(lubridate)
library(evd)
ann_max <- aggregate(TG ~  year(DATE),
                      FUN = max, data = data)
colnames(ann_max) <- c("Years", "TG")
```

J'ai récupéré les maxima annuels de températures à Agen et je fais l'inférence Bayesienne.

```{r}
nchain <- 10000
init <- fgev(ann_max$TG)$param
chain <- posterior(nchain, init, prior, lh = "gev",
                   data = ann_max$TG,
                   psd = c(1, 0.5, 0.1))
```

On affiche la chaine

```{r}
plot(mcmc(chain))
```

On peut alors fournir des estimations ponctuelles Bayesiennes pour les paramètres de la GEV

```{r}
summary(mcmc(chain))
```
