SIMULATION DE MONTE CARLO

Auteur·rice

TRA BI

Introduction

On considère le fleuve Coulantine traversant la ville d’Aussec. On souhaite savoir quelle est la probabilité que ce fleuve déborde et inonde la ville au cours de l’année 2024, selon qu’on construise ou non des murets autour de son lit. On modélise les pluies par un ensemble de dates, et des intensités pour chacune de ces dates. On considère ces pluies instantanées.

Notre modèle nous permetra de calculer la hauteur d’eau dans le fleuve à tout moment de l’année, et d’approximer, par la méthode de Monte Carlo, la probabilité que H_{max} la hauteur maximale de l’eau dans le fleuve dépasse le seuil h_0 (donné) correspondant à l’inondation de la ville.

Le présent rapport, récapitule la suite de fonctions et de techniques de programmation qui nous a permis d’estimer précisément P_0

Simulation des intensités

La variable aléatoire (I_n), l’intensité de la n\text{-ième} pluie, suit la loi de \mathbb{P}_I de densité f_I

f_I(x)=\mathbb{1}_{\mathbb{R}_{+}}(x)\frac{C_I}{m + \lvert x-x_0 \rvert^\lambda}

Les variables aléatoires (I_n)_{n\ge0} sont indépendantes identiquement distribuées (iid) et sont toutes d’une densité proportionelle à f(x)=\frac{1}{m + \lvert x-x_0 \rvert^\eta}

Code 
f <- function(x, p = params_def) {
  return(1/(p$m + abs(x - p$x0)^p$eta)) 
  }

cherchons une fonction de densité g telle que \frac{f(x)}{g(x)}\le c dont on connait la densité

remarque 🤓

f ressemble à la densité de Cauchy \quad \mathscr{C}(x_0,m) g(x)= \frac{m}{\pi}\frac{1}{m^2+(x-x_0)^2}

Code 
g <- function(x, p = params_def) {
  return(m/(pi*(p$m^2 + abs(x - p$x0)^2)))
  }

visualisons l’allure de la courbe de la fonction f

remarque 😲

la courbe de f a la forme d’une cloche, ce qui nous fait penser à la loi normal

Code 
f_sur_g <- function(x, p = params_def) {
  return((pi*(x - x0)^2)/(m^2 + m*abs(x - x0)^eta)) 
  }

visualisation

la fonction rintensite doit renvoyer un vecteur de taille n contenant les valeurs simulées suivant la loi \mathbb{P}_i

Simulons \mathbb{P}_i avec la methode de rejet en utilisant la loi Normal et Cauchy afin de choisir l’algorithme le plus éficace

  • Pour tester la validité de notre technique de simulation, on représente l’histogramme des réalisations obtenues par simulation surplantée de la courbe de la densité f
Code 
rintensite_Cauchy <- function(n, p = params_def){
   x <- numeric(n)
   afaire <- rep(TRUE, n)
   while(any(afaire)) {
     x[afaire] <- rcauchy(sum(afaire), p$x0, p$m )
     u <- runif(sum(afaire), min=0, max=10)
     afaire[afaire] <- f(x[afaire]) < u*(p$eta + 0.8)*dcauchy(x[afaire], location = p$x0, scale= p$m )
   }
   x
}

Code 
rintensite_Normal <- function(n, p = params_def){
   X <- numeric(n)
   afaire <- rep(TRUE, n)
   c <- p$eta + 1
   while(any(afaire)) {
     X[afaire] <- rnorm(sum(afaire), p$x0, p$m^1.5 )
     Y <- runif(sum(afaire), min=0, max=c)
     afaire[afaire] <- f(X[afaire]) < Y*dnorm(X[afaire], p$x0, p$m^1.5 )
   }
   X
}

Avec une bonne valeur de C_I, on observe une concordance entre la courbe de la densité f et les histogrammes obtenues. Cela nous rassure quant à la validité de notre technique de simulation.

durées de simulation de chaque méthode

Code 
benchmark(rintensite_Normal(10^3), rintensite_Cauchy(10^3))[,c(1,3,4)]
                     test elapsed relative
2 rintensite_Cauchy(10^3)    0.54     10.8
1 rintensite_Normal(10^3)    0.05      1.0

la simulation de la loi \mathbb{P}_i à partir de la loi normal est la plus éfficace.

Dans la suite rintensite_Normal sera rintensite

Code 
rintensite = rintensite_Normal

silmulaltion de p*

On modélise les dates des pluies par un processus de Poisson inhomogène, de fonction intensité \lambda. \lambda(t)=\lambda_0(1+\alpha\sin(4\pi t))

Code 
intens <- function(t, p = params){
  return(p$lambda0*(1 + p$alpha*sin(4*pi*t)))
  
}

\lambda_0(1+\alpha\sin(4\pi t)) \le 2\lambda_0 \quad \text{car } \begin{cases} 0 \le sin(4\pi t) \le 1 \\0 \le \alpha \le 1 \end{cases}

Soit 2\lambda_0 un majorant de la fonction λ

simulation des dates

on sait que :

  • 2\lambda_0 est un majorant de la fonction \lambda.

  • le temps t est le temps en années, avec t = 0 au début de l’année 2024 et t = 1 à la fin de l’année 2024.

étapes de la simulation:

  • Donc homo un processus de Poisson homogène d’intensité 2\lambda_0.

  • (U_i)_i des variables aléatoires indépendantes et indépendantes de (homo_i)_i, de loi U([0, 1]).

  • le vecteur d’elément croissant obtenu en triant les élémants du processus homogène homoh qui vérifie 2λ_0∗U_i \le λ(homoi) est un processus de Poisson inhomogène d’intensité λ.

Code 
Dates <- function(t = 1, p = params){
  N <- rpois(1, lambda = 2*p$lambda0*t)
  homo <- runif(N, min = 0, max= 1)
  U <- runif(N, min = 0, max = 1)
  inhomo <- sort(homo[2*p$lambda0*U <= intens(homo)])
  return(inhomo)
}

Les zones avec la pente plus élévée désignent les période où il pleut moins

Code 
a <- Dates()

barplot(a, col = rainbow(10))

r est la fonction de résorption de l’eau de ce fleuve r(t)=\mathbb{1}_{\mathbb{R}_{+}}(t)\exp(-vt)

Code 
r <- function(t){(t > 0)*exp(-v*t)}

La hauteur d’eau H(t) dans le fleuve Coulantine au temps t

H(t)=\displaystyle\sum_{i\ge1}^{}I_i\times\ r(t-T_i) les intensités suplémentaires qui n’ont pas de correspondance en dates ne concerne pas les pluie de l’année 2024.

nous n’en tiendront pas compte.

loi de poissson [0, 0.5] considère uniquement les pluies avant 0.5

Code 
Hauteur <- function(t){
  d <- Dates(t)
  H <- sum(rintensite_Normal(length(d))*r(t - d))
  return(H)
}

visualisation de l’évolution de la hauteur de l’eau du fleuve

[1] 189

on approxime 0\le t \le 1 par x = seq(from = 0, to = 1, by = 0.001)

ind.vector renvoie un vecteur binaire [0, 1], les 1 corespondent à des scénarios d’inondation

Code 
# à améliorer car la rapidité de petoile et de petoile.petite depend de lui

ind.vector <- function(n.simul, p = params){
   S <- rep(0, n.simul)
  t = seq(from = 0, to = 1, by = 0.01)
  for(i in 1:n.simul){
    x = 0
      if ( max(sapply(t, FUN = Hauteur) > p$h0) ){
        x = 1
      }
    S[i] = x
  }
  return(S)
}

probabilité d’inondation

on créé une liste de H_{max} et on détermine la propention de de ses éléments qui sont supérieur à h_0

Code 
petoile <- function(n.simul, p = params){
  S <- ind.vector(n.simul)
  probabilite = mean(S)
  return(list("p" = probabilite, "demi.largeur" = qnorm(0.95) * sd(S)/sqrt(n.simul)))
}
Code 
n.simul <- 10^3
system.time(a <- petoile(n.simul))
utilisateur     système      écoulé 
      37.13        0.55       39.62
Code 
a
$p
[1] 0.066

$demi.largeur
[1] 0.01292082

prob est la fonction qui renvoie l’évolution de l’estimation p en fonction du nombre de simulation

Code 
prob <- function(n.simul, p = params){
  S <- ind.vector(n.simul)
  return(cumsum(S)/1:n.simul)
}
Code 
n.simul = 1000
plot( 1:n.simul, y = prob(n.simul), t = 'l')

réduction de variance

I = E(h(X)) = P(\mathbb{1}_{H_{max}>h_0}) on a besoin d’une fonction l telle que E(l(X)) est facile à calculer Soit X la loi normal \mathscr{N(h_0, \sigma)}.

On choisit l(X) = 1_{X≥h_0} comme variable de contrôle.

On sait calculer son espérance facilement : E(l(X)) = E(1_{\mathscr{N(h_0, \sigma)} ≥ h_0})= \frac{1}{2}

On simule maintenant l’espérance de $ E(h(X) − l(X))$ par Monte-Carlo, pour obtenir une approximation de E(h(X)) = E(h(X) − l(X)) + E(l(X)).

(H.max(n.simul) >= seuil) = ind.vector|

Code 
petoile.petite <- function(n.simul, p = params){
  seuil <- p$h0
  x.norm <- rnorm(n.simul, mean = seuil, sd = 1)
  S <- ind.vector(n.simul) - (x.norm >= seuil)
  
  probabilite = mean(S) + 1/2
  return(list("p" = probabilite, "demi.largeur" = qnorm(0.95) * sd(S)/sqrt(n.simul)))
}
Code 
seuil <- p$h0
prob.petite <- function(n.simul, p = params){
  x.norm <- rnorm(n.simul, mean = seuil, sd = 1)
  S <- ind.vector(n.simul) - (x.norm > seuil)
  
  (cumsum(S) + 1/2)/1:n.simul
}

( en cours) tableau de comparaison des probas de petoile et petoile.petite avec des valeurs de h0 de plus en plus garndes

Utilisation d’un code C

(en cours) vu le TD qu’on a fait aujourd’hui, que chacun fasse la verson C++ des code ici present pour le vendredi.

ça me fera vraiment avancer le projet

Influence des paramètres

Code 
somme <- 0
for(i in 1:n.simul){
  somme <- somme + length(Dates())
}
somme/n.simul
[1] 170.275
  • \lambda_0 Si les dates de pluie suivait un processus de Poisson homogène d’intensité, \lambda_0. Le nombre moyen de pluie au cours l’année 2024 serait \lambda_0.

Mais nous somme en présence d’un processus de Poissons inhomogène d’intensité \lambda_0(1+\alpha\sin(4\pi t)). (1+\alpha\sin(4\pi t)) qu’on peut faire auxiller en fonction du parametre \alpha.

le nombre de pluie simulé par notre processus de poissons en 2024 est en moyenne \lambda_0 avec une fluctuation plus ou moins grande de \alpha

on calcule la moyenne du processus de poissons pour savoir il pleut environ combien de fois par an

  • lambda0 augmente le risque d’inondation
  • La distribution f des intensités est symétrique par rapport à x_0. x_0 est donc notre paramètre de position. On peut donc les considérer l’intensité moyenne des pluis.

  • x0 augmente l’intensité des pluies et donc fait croître rapidement le niveau de l’eau.

  • x0 augmente le risque d’inondation (futur graphe)

pour une valeur de v élevée, le niveau de l’eau du le fleuve baisse drastiquement.

(futur graphe)

v diminue donc le risque d’inondation

  • \alpha fait fluctuer la densité des pluies. il pleut plus pendant certaines période de l’année à l’instar des saisons pluivieuse.

plus \alpha est grand, plus les saisons pluvieuse sont plus longues.

  • \alpha augmente le risque d’inondation

  • TRIVIAL

  • le risque diminue voyons ! c’est évident !

  • ( on fera une courbe des proba en fonction de h0)

  • j’ai envie de laisser cette phrase pour faire rigoler les profs 🤣🤣