---
title: "Coagulation"
author: "Nous"
format: html
editor: visual
---

## Le code en R car c'est mieux

```{r}
coagulation <- function(nchain, data, init, prior){
  ## Création des objets de mes données
  time <- data$time
  diet <- data$diet
  ndiet <- length(unique(diet))
  freq <- table(diet)
  nobs <- length(time)

  ## Initialization
  chain <- matrix(NA, nchain, 3)
  current <- init
  names(current) <- colnames(chain) <- c("mu", "sigma2", "sigma2t")
  
  latent_chain <- matrix(NA, nchain, ndiet)
  latent <- rep(NA, ndiet)
  names(latent) <- colnames(latent_chain) <- unique(diet)
  
  ## Compute some fixed qunatitites once for all
  sum_diet <- aggregate(time ~ diet, data = data, FUN = sum)$time
  
  for (iter in 1:nchain){
    ## Update latent variables, i.e., theta_j
    cond_vars <- 1 / (1 / current["sigma2t"] + freq / current["sigma2"])
    cond_means <- cond_vars * (current["mu"] / current["sigma2t"] + sum_diet / current["sigma2"])
    latent[1:4] <- rnorm(ndiet, cond_means, sqrt(cond_vars))
    
    ## Update mu
    cond_var <- 1 / (1 / prior$mu[2] + ndiet / current["sigma2t"])
    cond_mean <- cond_var * (prior$mu[1] / prior$mu[2] + sum(latent) / current["sigma2t"])
    current["mu"] <- rnorm(1, cond_mean, sqrt(cond_var))
    
    ## Update sigma2
    cond_shape <- prior$sigma2[1] + 0.5 * nobs 
    cond_scale <- prior$sigma2[2] + 0.5 * sum((time - latent[diet])^2)
    current["sigma2"] <- 1 / rgamma(1, cond_shape, rate = cond_scale)
    
    ## Update sigma2t
    cond_shape <- prior$sigma2t[1] + 0.5 * ndiet
    cond_scale <- prior$sigma2t[2] + 0.5 * sum((latent - current["mu"])^2)
    current["sigma2t"] <- 1 / rgamma(1, cond_shape, rate = cond_scale)
    
    ## Update the chains
    chain[iter,] <- current
    latent_chain[iter,] <- latent
  }
  
  return(list(chain = chain, latent = latent_chain))
}
```

```{r}
A <- c(62, 60, 63, 59)
B <- c(63, 67, 71, 64, 65, 66)
C <- c(68, 66, 71, 67, 68, 68)
D <- c(56, 62, 60, 61, 63, 64, 63, 59)
data <- data.frame(time = c(A, B, C, D), diet = c(rep("A", length(A)), rep("B", length(B)), 
                   rep("C", length(C)), rep("D", length(D))))
```

```{r}
nchain <- 10^4
init <- c(mean(data$time), var(data$time), 1)
prior <- list(mu = c(0, 100), sigma2 = c(.1, 1), sigma2t = c(.1, 1))
ans <- coagulation(nchain, data, init, prior)
```

```{r}
latent_post_med <- colMeans(ans$latent)###apply(ans$latent, 2, median)
post_med <- colMeans(ans$chain)##apply(ans$chain, 2, median)
cols<- RColorBrewer::brewer.pal(4, "Set2")
i <- 1
plot(0, xlab = "Coagulation time", ylab = "Density",
     type = "n", xlim = c(50, 80), ylim = c(0, 0.2))
for (diet in unique(data$diet)){
  pdf <- function(x)
    dnorm(x, latent_post_med[diet], sqrt(post_med["sigma2"]))
  
  curve(pdf, from = 50, to = 80, col = cols[i], add = TRUE, lwd = 2)
  rug(data$time[data$diet == diet], col = cols[i], lwd = 2)
  i <- i + 1
}

pdf <- function(x)
  dnorm(x, post_med["mu"], sqrt(post_med["sigma2t"]))
curve(pdf, from = 50, to = 80, lwd = 2, add = TRUE)
legend("topright", c(unique(data$diet), "Sur-pop"), col = c(cols, "black"),
       bty = "n", lty = 1, lwd = 2)
```
