library(Rcpp)

Embedding Rcpp code in your R code

You can also write Rcpp code in your R code in 3 ways using sourceCpp() cppFunction() evalCpp() respectively.

sourceCpp()

Save Rcpp code as string object in R and compile it with sourceCpp()

src <-

"#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]

double rcpp_sum(NumericVector v){

  double sum = 0;

  for(int i=0; i<v.length(); ++i){

    sum += v[i];

  }

  return(sum);

}"



sourceCpp(code = src)

rcpp_sum(1:10)
[1] 55

cppFunction()

The cppFunction() offers a handy way to create a single Rcpp function. You can omit #include <Rcpp.h> and using namespase Rcpp; when you use cppFunction().

src <-

  "double rcpp_sum(NumericVector v){

    double sum = 0;

    for(int i=0; i<v.length(); ++i){

      sum += v[i];

    }

    return(sum);

  }

  "

Rcpp::cppFunction(src)

rcpp_sum(1:10)
[1] 55

evalCpp()

You can evaluate a single C++ statement by using evalCpp().

# Showing maximum value of double.

evalCpp('std::numeric_limits<double>::max()')
[1] 1.797693e+308

The table below presents the correspondence of data types between R/Rcpp/C++.

Value | R vector | Rcpp vector | Rcpp matrix | Rcpp scalar | C++ scalar |

|————|————|————-|————-|————|————|

Logical | logical | LogicalVector | LogicalMatrix | - | bool |

Integer | integer | IntegerVector | IntegerMatrix | - | int |

Real | numeric | NumericVector | NumericMatrix | - | double |

Complex | complex | ComplexVector | ComplexMatrix | Rcomplex | complex |

String | character | CharacterVector (StringVector) | CharacterMatrix (StringMatrix) | String | string |

Date | Date | DateVector | - | Date | - |

Datetime | POSIXct | DatetimeVector | - | Datetime | time_t |

data.frame, list

R | Rcpp |

|————–|————-|

data.frame | DataFrame |

list | List |

S3 class | List |

S4 class | S4 |

Warning: le package 'plotly' a été compilé avec la version R 4.3.2
Le chargement a nécessité le package : ggplot2
Warning: le package 'ggplot2' a été compilé avec la version R 4.3.2

Attachement du package : 'plotly'
L'objet suivant est masqué depuis 'package:ggplot2':

    last_plot
L'objet suivant est masqué depuis 'package:stats':

    filter
L'objet suivant est masqué depuis 'package:graphics':

    layout
Warning: le package 'rbenchmark' a été compilé avec la version R 4.3.1
Warning: le package 'shiny' a été compilé avec la version R 4.3.2
Warning: le package 'gapminder' a été compilé avec la version R 4.3.2
Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Warning: `line.width` does not currently support multiple values.

Introduction

Explanation: in a quarto doc, the ::: operator adds a html div. Then, writing .panel-tabset .nav-pills adds 2 classes to the div. Quarto relies on a framework called bootstrap. Bootstrap knows thoses classes and transform their children in tabs automatically for us!

Here is the link to the section 2 of the document!

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\)

Oganisation

  • Metushela : expert en chatGpt
  • La reine mère
  • Pachelle l’intello
  • Tra Bi : chargé des

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^\lambda}\)

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}\]

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

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

visualisation

x = seq(from = 0, to = 4*p$x0, by = 0.001)

plot(x, f(x), type = 'l', col= 'blue', main = "f")
lines(x, (p$eta + 0.8)*dcauchy(x, location = p$x0, scale= p$m ) , col = 'red', main = "cauchy")
lines(x, (p$eta+0.7 )*dnorm(x, mean = p$x0, sd = p$m^1.5 ), col = 'green', main= "exp")
legend("topright", legend=c("normale", "cauchy", "f"), col=c("green", "red", 'blue'), lty=c(1, 1))

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\)

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

benchmark(rintensite_Normal(10^3), rintensite_Cauchy(10^3))[,c(1,3,4)]
                     test elapsed relative
2 rintensite_Cauchy(10^3)    1.86     7.75
1 rintensite_Normal(10^3)    0.24     1.00

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

Dans la suite rintensite_Normal sera rintensite

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))\]

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é λ.

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

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)\]

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

Hauteur1 <- function(t, p = params){
  Ti <- Dates(t)
  H <- sum(rintensite_Normal(length(Ti))*r(t - Ti))
  return(H)
}

pourquoi pour chaque t, on referai la simulation ? on pourait fixer les dates de pluies à l’avance, biensur par la simulation de poisson pour respecter l’aléa. et pour chaque instant \(t\) considérer que les pluies qui sont déjà survenu.

Pour un instant de \(t\) donné, seuls les instants \([T_1,...,T_n]\) tel que \(T_n\le t<T_{n+1}\) seront considéré, alors il y aura eu n pluies avant cette date.

Hauteur <- function(t, p = params){
  Ti <- Dates(t, p)
1  Ti <- Ti[ Ti <= t ]
2  Ii <- rintensite_Normal(length(Ti) )
  H <- sum(Ii*r(t - Ti))
  return(H)
}
1
les pluies survenues après \(t\) n’intervienent pas dans le calcul de la hauteur.
2
On simul les intensité pour chaque date de pluie

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

nous allons visualiser l’évolution de la hauteur de l’eau sur 366 jours car 2024 est une année bixectille

# Libraries
library(ggplot2)
library(dplyr)

Attachement du package : 'dplyr'
Les objets suivants sont masqués depuis 'package:stats':

    filter, lag
Les objets suivants sont masqués depuis 'package:base':

    intersect, setdiff, setequal, union
t = seq(from = 0, to = 1, length.out = 366 )


# Dummy data
data <- data.frame(
  day = as.Date("2024-01-01") + 1:length(t),
  value = sapply(t, FUN = Hauteur1)
)

# Most basic bubble plot
p <- ggplot(data, aes(x=day, y=value)) +
  geom_line(color="#69b3a2") +
  geom_area(fill = "#69b3a2", alpha = 0.5)+
  xlab("Date")+
  ylab("niveau de l'eau") # +  theme_ipsum()

ggplotly(p)

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

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

ind.vector1 <- 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)
}
ind.vector <- function(n.simul, p = params){
   S <- rep(0, n.simul)
  t = c(seq(from = 0.6 , to = 0.9, by = 0.03) )
  for(i in 1:n.simul){
    x = 0
      if ( max(sapply(t, FUN = Hauteur) > p$h0) ){
        x = 1
      }
    S[i] = x
  }
  return(S)
}
ind.vector1 <- function(n.simul, p = params){
  t <- seq(from = 0, to = 1, by = 0.01)
  data <- rep(t, n.simul)
  M <- matrix(sapply(data, FUN = Hauteur), ncol = n.simul)
  S <- apply(M, 2, max) > p$h0
  return(S)
}
system.time(ind.vector(10^3))
utilisateur     système      écoulé 
       3.15        0.20        7.97 

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\)

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)))
}
n.simul <- 10^3
system.time(a <- petoile(n.simul))
utilisateur     système      écoulé 
       3.10        0.06        7.66 
a
$p
[1] 0.037

$demi.largeur
[1] 0.009823325

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

prob <- function(n.simul, p = params){
  S <- ind.vector(n.simul)
  return(cumsum(S)/1:n.simul)
}
# Get data:
library(gapminder)
 
# Charge libraries:
library(ggplot2)
library(gganimate)
 
# Make a ggplot, but add frame=year: one image per year
ggplot(gapminder, aes(gdpPercap, lifeExp, size = pop, color = continent)) +
  geom_point() +
  scale_x_log10() +
  theme_bw() +
  # gganimate specific bits:
  labs(title = 'Year: {frame_time}', x = 'GDP per capita', y = 'life expectancy') +
  transition_time(year) +
  ease_aes('linear')

# Save at gif:
#anim_save("271-ggplot2-animated-gif-chart-with-gganimate1.gif")
n.simul = 10000
df <- data.frame(a = 1:n.simul, b = prob(n.simul) )
ggplot(df, aes(x = a, y = b))+
  geom_line()

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|

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)))
}
#petoile.petite(1000)
seuil <- p$h0
prob.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)
  
  (cumsum(S)/(1:n.simul) + 1/2)
}

( 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

Quelles sont les techniques de programmation qui vous semblent intéressantes dans votre code ?

processus de poissons methode de rejet 2 reduction de variance

Quelle est la partie du code que vous avez codée en C/C++, et pourquoi ?

c_intensité c_Dates c_hauteur c_ind.vecteur

Comment avez-vous vérifié la validité de votre code (tests avec certaines valeurs de paramètres, calculs intermédiaires, représentations graphiques, etc.) ?

Quelles études basées sur ce modèle vous paraîtraient pertinentes ? Quelles évolutions de ce modèle proposeriez-vous, afin de le rendre plus réaliste ?