Monte Carlo :

Bertrand PATUREL & Omar El Mellouki

Professeur Pierre Gloaguen

Cours ENSTA : https://papayoun.github.io/courses/2020_monte_carlo/

Dans ce problème on souhaite savoir si à partir d'un vecteur d'observations $X_{1:n}=(X_1,..,X_n)$ il existe une suite d'observations (noté $S(i,j)=\sum_{k=i}^j{R_k}$) anormalement élevées. Nous testons donc l'hypothèse $H_0$ : le vecteur est i.i.d. contre $H_1$ il existe une fenêtre de valeur anormalement élevées.

Calcul des probabilités sous l'hypothèse $H_0$.

Soit $k$ un en entier naturel compris entre 1 et n. $R_k$ est le rang de $X_k$ cela signifie que $X_k$ est plus grand que $R_k-1$ variable(s) aléatoire(s) du vecteur $X_{1:n}$. Donc $R_k=\sum_1^n\mathbb{1} _{X_k \geq X_i}$. Or pour chaque $i$ différent de $k$, $\mathbb{1} _{X_k \geq X_i}$ suit une loi de bernouilli de paramètre $p=\frac{1}{2}$, puisque l'on a autant de chance d'obtenir $X_k \geq X_i$ que son évènement complémentaire. La variable $R_k-\mathbb{1} _{X_k \geq X_i}$ somme de $n-1$ variables indépendantes de Bernouilli de paramètre $p=\frac{1}{2}$ suit une loi binomiale $\mathbb{B}(n-1,\frac{1}{2})$, donc une loi uniforme sur {0,..,n-1}. Puisque l'on rajoute $1=\mathbb{1} _{X_k \geq X_i}$ à la variable $R_k-\mathbb{1} _{X_k \geq X_i}=\sum_{i=1, i \ne k}^n\mathbb{1} _{X_k \geq X_i}$, $R_k$ suit alors une loi uniforme sur {1,...,n}.

Soit $(i,j)$ couple d'entiers naturels telque $i\ne j$ compris entre 1 et n. Alors on souhaite calculer $\mathbb{P}(R_k=i | R_l=j)=\frac{\mathbb{P}(R_k=i\cap R_l=j)}{\mathbb{P}(R_l=j)}$ :

  • $\mathbb{P}(R_k=i\cap R_l=j)$ est en combinatoire "la probabilité de placer $X_k$ en i ème position (dans la ième case) ou il y a $n$ choix de positions (i.e. de cases), puis de placer $X_l$ en jème poisition (dans la jème case) ou il y a $n-1$ choix de positions restantes (i.e. de cases restantes)". On a donc une probabilité de $\mathbb{P}(R_k=i\cap R_l=j)=\frac{1}{n(n-1)}$.

  • $\mathbb{P}(R_l=j)$ est en combinatoire "la probabilité de placer $X_l$ en j ème position (dans la jème case) parmis les $n$ choix de positions (i.e. de cases)". On a donc une probabilité de $\mathbb{P}(R_l=j)=\frac{1}{n}$.

  • $\mathbb{P}(R_k=i | R_l=j)$ est en combinatoire "la probabilité de placer $X_k$ en i ème position (dans la ième case) sachant que $X_l$ est deja en jème position (dans la jème case)", il reste donc $n-1$ choix de positions (i.e. de cases) pour placer la variable $X_k$. On a donc une probabilité de $\mathbb{P}(R_k=i | R_l=j)=\frac{1}{n-1}$.

Donc $\mathbb{P}(R_k=i | R_l=j)=\frac{1}{n-1}$.

$S(i,j)$ est une fenêtre temporelle des rangs des variables comprises entre $X_i$ et $X_j$ incluses. Cette variable aléatoire est naturellement définit par la somme des rangs inclus dans cette fenêtre. $S(i,j)$ va nous permettre de détecter les plages de valeurs anormalement élevées (si elles existent) du vecteur $X_{1:n}$. Pour ce faire nous allons la comparer à sa valeur de référence (sa moyenne) et étudier son écart type (en calculant sa variance).

Soit $(i,j)$ un couple d'entier entre 1 et n.

Par définition : $\mathbb{E}[S(i,j)]=\mathbb{E}[\sum_{k=i}^j{R_k}]$

Par linéarité : $\mathbb{E}[S(i,j)]=\sum_{k=i}^j\mathbb{E}[{R_k}] = (j-i+1)\mathbb{E}[R_1] = \frac{(j-i+1)(n+1)}{2}$ puisque l'espérance de la variable aléatoire $R_1$ est celle de la loi uniforme sur {1,..,n} : $\mathbb{E}[R_k]=\mathbb{E}[R_1]=\frac{n+1}{2}$

Soit $(i,j)$ un couple d'entier entre 1 et n.

Par définition : $\mathbb{V}[S(i,j)]=\mathbb{V}[\sum_{k=i}^j{R_k}]$

Commencons par calculer la variance de $R_1$ qui correspond à celle de toutes les variables $R_k$ suivant la même loi : $\mathbb{V}[R_1]=\mathbb{E}[R_1^2]-\mathbb{E}[R_1]^2=\sum_{l=1}^n l^2 \mathbb{P}(R_1=l)-\frac{(n+1)^2}{4}=\frac{\sum_{l=1}^n l^2}{n}-\frac{(n+1)^2}{4}=\frac{2n^2+3n+1}{6}-\frac{(n+1)^2}{4}=\frac{n^2-1}{12}=\frac{(n-1)(n+1)}{12}$

Nous pouvons aussi calculer la covariance de $R_1$ et $R_2$ qui est correspond à celle de tous les couple $R_k$ et $R_l$ : $\mathbb{Cov}(R_k,R_l)=\mathbb{Cov}(R_1,R_2)=\mathbb{E}[R_1R_2]-\mathbb{E}[R_1]\mathbb{E}[R_2]=\mathbb{E}[R_1R_2]-\frac{(n+1)^2}{4}=\sum_{1 \le p \ne l \le n}lp\mathbb{P}(R_1R_2=lp)-\frac{(n+1)^2}{4}=2\sum_{1 \le p < l \le n}lp\mathbb{P}(R_1=p\cap R_2=l)-\frac{(n+1)^2}{4}=2\sum_{1 \le p < l \le n}lp\frac{1}{n(n-1)}-\frac{(n+1)^2}{4}=2\sum_{1 \le p < l \le n}lp\frac{1}{n(n-1)}-\frac{(n+1)^2}{4}=2\frac{1}{n(n-1)}\sum_{1 \le p \le n}\sum_{p+1 \le l \le n}lp-\frac{(n+1)^2}{4}=2\frac{1}{n(n-1)}\sum_{1 \le p \le n}l\sum_{p+1 \le l \le n}p-\frac{(n+1)^2}{4}=2\frac{1}{n(n-1)}\sum_{1 \le p \le n}l(\frac{n(n+1)}{2}+\frac{l(l+1)}{2})-\frac{(n+1)^2}{4}=[2\frac{1}{n(n-1)}(\sum_{1 \le p \le n}\frac{-l^2}{2}+\frac{-(l+1)^2}{2})+(\frac{n(n+1)}{2})^2]-[\frac{(n+1)^2}{4}]=[\frac{3n^2+5n+2}{12}]-[\frac{3n^2+6n+3}{12}]=\frac{-n-1}{12}$

Par propriété de la variance : $\mathbb{V}[S(i,j)]=[\sum_{k=i}^j\mathbb{V}[{R_k}]]+[2\sum_{i \le p < l \le j}\mathbb{Cov}(R_p,R_l)]=[\sum_{k=i}^j\mathbb{V}[{R_k}]]+[2\sum_{i \le p < l \le j}\mathbb{Cov}(R_1,R_2)]=[(j-i+1)(\mathbb{V}[R_1])]+[(j-i+1)(j-i)(\mathbb{Cov}(R_1,R_2))]=[(j-i+1)(\frac{(n-1)(n+1)}{12})]+[(j-i+1)(j-i)(\frac{-n-1}{12})]=\frac{j-i+1}{12}((n+1)(n-1)+(n+1)(i-j))=\frac{j-i+1}{12}(n+1)(n+i-j-1)$

Choisir entre les deux hypothèses par Monte Carlo.

Pour simuler $T_n$ sous $H_0$, on va générer une réalisation de n variables aléatoires i.i.d. Leur loi n'est pas importante, tant qu'ils sont i.i.d, et on sait que la variable aléatoire $R_k$ va dans tout les cas suivre une loi uniforme sur $[1:n]$. On va ensuite calculer la matrice triangulaire des $S(i,j)$, puis celle des $T(i,j)$ à l'aide des expressions de $m_{ij}$ et $v_{ij}$, en faisant attention à bien remplacer la valeur de $T(1,n)$ par 0. On retrouve finalement la valeur maximale de la matrice $T$, qui va correspondre à la valeur simulée de notre statistique.

On souhaite déterminer le quantile $q_{1-\alpha}$ de la loi de $T_n(X)$ sous l'hypothèse $H_0$.

Ce quantile de niveau $\alpha$ vérifie : $\mathbb{P}(T_n(X)<q_{1-\alpha})=\alpha$.

Nous sommes de capable de calculer $\mathbb{P}(T_n(X)<q)=\mathbb{E}[\mathbb{1}_{T_n(X)<q)}]=\frac{1}{M}\sum_{l=1}^{M}\mathbb{1}_{T_n^l(X)<q}$ par méthode de Monte Carlo pour un q donné en générant M réalisations de $T_n$. Puis on on détermine, par une boucle for le quantile $q_{1-\alpha}$ compris entre 0 et 1 vérifiant $\mathbb{P}(T_n(X)<q_{1-\alpha})=\frac{1}{M}\sum_{l=1}^{M}\mathbb{1}_{T_n^l(X)<q_{1-\alpha}} \approx \alpha $.

On peut pour se faire déterminer $q_{1-\alpha}$ de la forme $0.01k$ ou $k$ est compris entre $0$ et $100$. On choisit $q_{1-\alpha}$ telque $q_{1-\alpha}= $argmax$ (0.01k | \mathbb{P}(T_n(X)<0.01k) \leq \alpha)$. Sinon on peut trouver grace à la fonction uniroot le zéro de la fonction $f(q) = \mathbb{P}(T_n(X)<q)-\alpha = \frac{1}{M}\sum_{l=1}^{M}\mathbb{1}_{T_n^l(X)<q} - \alpha$ ou encore la fonction quantile à partir du jeu de données.

Ensuite nous comparons le quantile trouvé à $t*$ notre statistique obtenue sur le jeu de donnée, si $t*$ est supérieur au quantile, alors on rejette l'hypotèse $H_0$.

Nous pouvons confirmer notre choix en calculant la p-value à l'aide de la méthode de Monte Carlo : On approxime par l'estimateur de Monte Carlo la p-value $\mathbb{P}(T_n(X)<t*=\mathbb{E}[\mathbb{1}_{T_n(X)<t*)}]\approx \frac{1}{M}\sum_{l=1}^{M}\mathbb{1}_{T_n^l(X)<t*}$ ou $t*$ est la statistique calculée sur le jeu de données des températures de Hobart. Si cette p-value est inférieur à $\alpha$ alors on rejette $H_0$.

Implémentation du code sur un jeu de données à tester.

In [14]:
hobart <- read.table("https://papayoun.github.io/courses/2020_monte_carlo/enonces_td/hobart.txt
",sep = ";", header = TRUE)# On importe le jeu de données des températures de Hobart 

hobart = hobart$mean_january_temperature

get_tn = function(X,a=0){
    
    n = length(X)
    
    #Le vecteur des Rk
    R = order(X,decreasing= T)
  
    #La matrice T(i,j)
    T_ = matrix(data = 0,nrow = n,ncol=n)
    for(j in 1:n){
      for(i in 1:j){
        if (i!=1 & j!=n){
        T_[i,j] = (sum(R[i:j]) - (n+1)*(j - i +1)/2)/sqrt(((1/12)*(n+1)*(j - i+1)*(n - j+i - 1))) 
          }
        }
    }
    
    #Le cas limite
    T_[1,n] = 0
    
    #Le résultat avec le max et les indices
    res = list(max(T_),which(T_ == max(T_), arr.ind = TRUE))
    
    #On renvoie la valeur maximale t* et si a vaut 1 on renvoie l'emplacement
    if(a == 1){return(res)}
    if(a == 0){return(max(T_))}
    
    
}

(t_star=get_tn(hobart)) #On affiche t*
3.93283106029876
In [15]:
get_ho_sample = function(n,M){
  res = c(1:M)
  for(a in 1:M){
    
    #Le vecteur des Xk
    X = rnorm(n,0,1)
    
    res[a] = get_tn(X,0)
  }
  return(res)
}
get_ho_sample(50,100) #On affiche notre sample de Tn pour 100 valeurs générées de Tn.
  1. 2.89655408640121
  2. 2.29697157826834
  3. 1.85869176721789
  4. 2.71672052482072
  5. 2.14464434678988
  6. 2.8062827399247
  7. 3.14907393794853
  8. 2.51229766338243
  9. 1.98029508595335
  10. 3.65154764267719
  11. 2.44960981286696
  12. 2.42399304803433
  13. 2.40049471568299
  14. 2.21711547354167
  15. 1.73614066168704
  16. 2.39591510850471
  17. 2.56689220152107
  18. 2.7649061254142
  19. 2.53853954204695
  20. 2.17832459454868
  21. 2.70022995873784
  22. 2.32924535991426
  23. 2.30226413950859
  24. 2.02980246310218
  25. 2.72566903872856
  26. 3.52190689796676
  27. 2.35259556285243
  28. 3.2540938154446
  29. 2.71775445143818
  30. 2.12881721739985
  31. 3.11415742546652
  32. 2.55825234707815
  33. 2.76894838583147
  34. 2.15856706282336
  35. 2.12881721739985
  36. 2.26719545232073
  37. 2.99127270878144
  38. 2.43059692636186
  39. 2.54021365676096
  40. 2.46743496248992
  41. 2.02209324125903
  42. 2.31192661400272
  43. 2.78388721237131
  44. 2.75351437843078
  45. 3.01526124183399
  46. 1.73275820020918
  47. 2.52118601360599
  48. 2.37707525504218
  49. 2.17832459454868
  50. 2.21711547354167
  51. 2.12881721739985
  52. 2.59399840040299
  53. 2.86599088332735
  54. 2.46205056200535
  55. 2.28863532752688
  56. 2.43167503549732
  57. 2.05892856874907
  58. 2.02980246310218
  59. 2.33458001072159
  60. 2.30681687311975
  61. 2.93231401339382
  62. 2.86607937218059
  63. 2.44152529203242
  64. 2.91691675598474
  65. 2.65700508925572
  66. 2.56745016631409
  67. 2.14464434678988
  68. 2.44803387950878
  69. 2.52142521985461
  70. 2.02209324125903
  71. 1.98029508595335
  72. 2.00751815055821
  73. 3.08977108226532
  74. 2.71775445143818
  75. 2.47144729487214
  76. 3.163101296733
  77. 2.43167503549732
  78. 2.182820625327
  79. 2.22634508381044
  80. 2.72600659934557
  81. 2.86079415940861
  82. 2.05346829197429
  83. 1.94039250423846
  84. 2.63555379206149
  85. 2.4738633753706
  86. 2.47144729487214
  87. 2.07930984025102
  88. 2.79825024295441
  89. 2.22783197169752
  90. 2.14860386108807
  91. 2.44152529203242
  92. 2.37674698635794
  93. 2.2555813128379
  94. 2.40918720869424
  95. 2.17037412052756
  96. 2.22634508381044
  97. 1.83177295450685
  98. 2.22783197169752
  99. 2.38974655785158
  100. 2.8062827399247

Calcul de l'estimation de la p value

In [20]:
set.seed(123)
n = length(hobart)

sample_hobart = get_ho_sample(n,5000)

#Fonction indicatrice(Tn>t*)

indicatrice = function(x){
  value = 0
  if(x>t_star){value = 1}
  return(value)
}

#Monte Carlo pour le calcul de E(indicatrice(Tn>t*))

indicatrice_Tn = lapply(sample_hobart,indicatrice)


I_hat <- mean(unlist(indicatrice_Tn)) # Estimateur de Monte Carlo de la proba 
gamma2_hat <- var(unlist(indicatrice_Tn)) # Variance de Monte Carlo de la proba
z_975 <- qnorm(0.975, 0, 1) # Quantile à 0.975 de la loi normale centrée réduite
IC = I_hat + z_975 * sqrt(gamma2_hat / 5000) * c(-1, 1) #Intervalle de confiance


### On affiche finalement l'intervalle de confiance à 95 %
library(ggplot2)
library(dplyr)

alpha=0.05 #risque de première espèce
my_samples = tibble(index = 1:5000, sample = unlist(indicatrice_Tn))
my_samples %>% 
  mutate(Estimate = cumsum(sample) / index, 
         IC_95_inf = Estimate - z_975 * sqrt(gamma2_hat / index), 
         IC_95_sup = Estimate + z_975 * sqrt(gamma2_hat / index)) %>%
  ggplot(aes(x = index)) +
  geom_ribbon(aes(ymin = IC_95_inf, ymax = IC_95_sup), 
              fill = "lightblue") + 
  geom_line(aes(y = Estimate)) + 
  geom_line(aes(y = alpha),color='red')#le risque de première espèce est en rouge
  labs(x = "M", title = "Estimation Monte Carlo de I",
       y = expression(hat(I)[M])) + 
  theme_bw() 

I_hat#Voici la p-value approchée, elle est inférieur au seuil alpha

NULL
0.0088

Q10

In [17]:
#On génère un vecteur contenant 5000 réalisations de Tn sous H0 de longueur n égale à celle de Hobart
vect = get_ho_sample(100,5000)  
#On calcule son quantile ) 0.95
quantile(vect,0.95) # Voici le quantile à 95%
t_star #Voici la statistique t*
95%: 3.49721352060894
3.93283106029876

On voit donc bien que le quantile est plus faible que la statistique $t*$ observée sur le jeu de données. On rejette donc $H_0$. Nous pouvions déjà rejeter $H_0$ en question 9 puisque la p-value (égale à environ 0.0088) est plus faible que $\alpha = 0.05$. Ainsi le jeu de données contient une fênetre de valeurs anormalement élevées.

Remerciements

Ce problème nous a été donné par Mr Pierre Gloaguen à l'ENSTA Paris : https://papayoun.github.io/courses/2020_monte_carlo/ .