1.1 Conseils pour les débudants de R

Si vous n'avez jamais suivi un cours de R, vous trouvrez dans les liens suivants deux documents par madame Tabea Rebafka pour le cours 4M015.

• Prise en main de R
• Programmation en R

2.1 Matrices des figures

Cas simple :

par(mfrow=c(1,3))
hist(rnorm(500), breaks=30)
hist(rexp(500), breaks=30)
hist(runif(500), breaks=30)

Cas (un peu) compliqué:

# One figure in row 1 and two figures in row 2
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
hist(rnorm(500), breaks=30)
hist(rexp(500), breaks=30)
hist(runif(500), breaks=30)

2.2 lines et curve

hist(rexp(300), freq=FALSE, breaks = 40)
x<-seq(0,8,by=.1)
# ici, x est un vecteur
lines(x, dexp(x), col = 'darkred')

hist(rexp(300), freq=FALSE, breaks = 40)
# ici, x est qqch abstrait
# pour aller plus loin  ?expression
curve(dexp(x), col = 'darkred', add=TRUE, from=.001)

2.3 list to data.frame

a <- list()
b <- list() 
a$X <- 1
a$Y <- 2
b$X <- 3
b$Y <- 4
data.frame(do.call(rbind, list(a,b)))

##   X Y
## 1 1 2
## 2 3 4

2.4 Visualisation 3D par plotly et ggplot2

library(ggplot2)

a <- seq(-10,10,0.5)
b <- a

f <- function(x,y){
      return(exp(-.01 *((x-1)^2 + (x+1)^2 + y^2)))
}

# z est une matrice t.q. z(i,j) = f(a_i,b_j) où
# (a_i,b_j)est le maillage. 

### par plotly
# library(plotly)
# plot_ly(x = a, y = b, z = sapply(1:length(b), function (i) return(f(a,b[i]))), type = 'surface')


### par ggplot2
A = length(a)
B = length(b)
df <- data.frame(a=rep(a,times=B),b=rep(b,each=A),f_z=c(sapply(1:length(b), function (i) return(f(a,b[i])))))
ggplot(df, aes(a, b, z=f_z)) + geom_contour(bins = 60, aes(colour = ..level..))

3.1 un exemple de optimize

On remarque que cette fonction n’est pas (du tout) magique pour l’optimisation non-convexe. Mais dans le cas convexe, elle marche très bien.

f<-function(x){return(x^2+sin(5*x))}
opt<-optimize(f, interval=c(-3,3), maximum=F, tol=0.01)
x_min<-opt$objective

curve(f(x), from=-3, to=3)
points(x_min, f(x_min), col='darkred',pch=19)
grid()

3.2 un exemple de findInterval

rappel de la définition de l’inverse généralisée d’une fonction réelle \(F\), noté \(F^{-1}\)

\[ F^{-1}(x) = \inf\{y\in \mathbb{R}: F(y)\geq x\} \]

x_inv<-findInterval(0.6, cumsum(dbinom(0:10,10,0.2)))

plot(0:10,cumsum(dbinom(0:10,10,0.2)))
points(x_inv,cumsum(dbinom(0:10,10,0.2))[x_inv+1] , col='darkred',pch=19)
grid()

3.3 Copule Gaussienne

On commence par un cas particulier (\(dim = 2\)). La copule est un objet pour décrire la dépendance entre les v.a. Pour les v.a. gaussiennes, la dépendance est totalement définie par la matrice de covariance.

Soit \((X_1,X_2) \sim \mathcal{N}_2 (0, \Gamma_2)\), où la matrice de covariance est

\[\Gamma_2 = \begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix}\]

Par construction, on a bien

\[\mathcal{L}(X_1) = \mathcal{L}(X_2) = \mathcal{N}(0,1).\]

Alors, comme la f.d.r. \(\Phi\in \mathscr{C}^0\)

\[(U_1,U_2) = (\Phi(X_1),\Phi(X_2))\]

est la copule qu'on cherche à simuler.

On revient au cas où \(dim = n\). Si on simule \((Z_1,\dots,Z_n)\) par

\[\forall i \in \{1,2,\dots,n\}\qquad Z_i = \sqrt{\rho} G_0 + \sqrt{1-\rho} G_i\]

où \(G_i\) sont les normales standards indépendantes, il est facile de vérifier que

\[\mathrm{Var}[Z_i] = 1\qquad et \qquad \mathrm{Cov}[Z_i Z_j] = \rho \qquad \forall i\neq j\in\{1,2,\dots,n\}\]

Alors, la matrice de covariance associée est de la forme

\[ \Gamma_n = \begin{pmatrix} 1 & \rho & \cdots & \rho \\ \rho & 1 & \cdots & \rho \\ \vdots & \vdots & \ddots & \vdots \\ \rho & \rho & \cdots & 1 \end{pmatrix} \]

Donc, pour simuler \((U_1,U_2,\dots,U_n)\) la copule gaussienne associée, il faut juste prendre

\[U_i = \Phi(Z_i) \qquad \forall i \in \{1,2,\dots,n\}\]

L’intérêt de cette méthode est que si vous voulez simuler les vecteurs gaussiens avec \(\Gamma\) différents, on peut simuler \((G_i)_{1\leq i \leq N}\) une seule fois. Ceci nous permet aussi de faire la simulation ‘dépendante’ parallèlement.

3.4 Urne d'Ehrenfest

On remarque que la loi invariante est bien \(\mathcal{B}(d,\frac12)\).

\[\forall k \in [0:d], \qquad \mathbf{P}(X_n = k) = \frac{1}{2^d}\binom{d}{k}\]

3.5 Processus de Poisson inhomogène

Soit \(\lambda:[0,+\infty[ \to ]0, +\infty[\) une fonction strictement positive, appelée fonction d’intensité. On considère maintenant des instants de sauts \((T_n)_{n \geq 1}\) (avec toujours \(T_0 = 0\)) donnés par la dynamique \(T_{n} = T_{n-1} + S_n\) où pour tout \(n \geq 1\), \(S_n\) est une variable aléatoire sur \(]0,+\infty[\) définie conditionnellement à \(T_{n-1}\) par

\[ \begin{aligned} \forall s > 0, \forall t_{n-1} > 0, \quad \mathbf{P} \bigl[S_n > s \; \bigl\vert \; T_{n-1} = t_{n-1} \bigr] &= \exp \bigl(-\int_0^s \lambda(u + t_{n-1}) d u \bigr) \\ &= \exp \bigl(- (\Lambda(t_{n-1} + s) - \Lambda(t_{n-1})) \bigr) \\ \end{aligned} \]

avec \(\Lambda(t) = \int_0^t \lambda(u) d u\) l’intensité intégrée. Dans la suite on note \(\Lambda^{-1}\) la réciproque de \(\Lambda\).

3.5.1 Simulation par inverse de l'intensité intégrée

Montrer que si \(E \sim \mathcal{E}(1)\) et que \(t_{n-1}\) est fixé, alors la variable aléatoire \(X = \Lambda^{-1}\bigl(E + \Lambda(t_{n-1})\bigr) - t_{n-1}\) vérifie

\[ \forall s > 0, \quad \mathbf{P} \bigl[X > s \bigr] = \exp \bigl(- (\Lambda(t_{n-1} + s) - \Lambda(t_{n-1})) \bigr), \]

Preuve :

Puisque \(\lambda\) est une fonction strictement positive, l’application \(\Lambda : t \rightarrow \Lambda(t)\) est continue et strictement croissante sur \(\mathbb{R}_+\) (donc bijective).

Alors, on a

\[ \begin{aligned} \mathbf{P}\big[ X > s\big] &= \mathbf{P}\big[\Lambda^{-1}\big(E + \Lambda(t_{n-1})\big) - t_{n-1} > s \big]\\ &=\mathbf{P}\big[E + \Lambda(t_{n-1}) > \Lambda(t_{n-1}+s)\big]\\ &=\mathbf{P}\big[E > \Lambda(t_{n-1}+s)-\Lambda(t_{n-1}) \big]\\ &=\exp \bigl(- (\Lambda(t_{n-1} + s) - \Lambda(t_{n-1})) \bigr) \end{aligned} \]

C’est-à-dire,

\[\mathcal{L}\big(X \; \big\vert \; T_{n-1} = t_{n-1}\big) = \mathcal{L}\big(S_n \; \big\vert \; T_{n-1} = t_{n-1}\big)\]

□
En déduire une méthode de simulation des instants de sauts \((T_n)_{n \geq 1}\) et du processus de comptage associé (appelé processus de Poisson inhomogène de fonction d’intensité \(\lambda\)).

Tester cette méthode avec la fonction d’intensité \(\lambda(t) = 0.1 + 5t\).

lambda <- function(t) return(0.1 + 5*t)
Lambda <- function(t) return(0.1*t + 2.5*t^2 )
Lambda_inv <- function(s) return(-0.02 + 0.2*sqrt(0.01+10*s))

X_cond <- function(t) {
    E <- rexp(1);
    X <- Lambda_inv (E + Lambda(t)) - t;
    return(X);
}


# svglite("poisson_inhomogene.svg")
Npaths <- 200
T <- 4
plot(NULL, xlim=c(0,T), ylim=c(0,25), xlab="Time t", ylab="Processus de comptage N_t",
     main = paste(Npaths, "trajectoires d'un processus de Poisson inhomogène"))

# on fait simple et on illustre ici l'utilisation d'une boucle for et d'une boucle while
for (k in 1:Npaths) {
    T_n <- 0
    path <- 0
    while (T_n < T) {
        T_n <- T_n + X_cond(T_n) 
        if (T_n < T) path <- c(path, T_n)
    }
    N_T <- length(path)-1
    N <- 0:N_T

    # on ajoute pour faire joli le dernier point (T, N_T)    
    path <- c(path, T)
    N <- c(N, N_T)
    
    lines(path, N, type="s", col=rgb(0,0,1,alpha=0.5))
}
grid()

tau_update <- function(tau,lambda_bar, T_bar){
    
    # tau: tau_n
    # lambda_bar: la borne de l'indensité
    # T_bar: un vecteur des instants de sauts d'un processus de Poisson
    
    k <- tau
    while(runif(1) > T_bar[k]/lambda_bar){
        k <- k+1
    }
    # la sortie: tau_{n+1}
    return(k)
}

3.6 Variable de contrôle

Soit \(X\) une v.a. et on veut estimer sa moyenne \(\mathbf{E}[X]\) par la méthode de Monte-Carlo. Soit \(Y\) une v.a. telle que l’on connaît la moyenne \(\mathbf{E}[Y]\). Alors, on a pour, tout \(\alpha\),

\[\mathbf{E}[X] = \mathbf{E}[X-\alpha Y] + \alpha \mathbf{E}[Y]\]

On suppose que l’on sait simuler suivant la loi \(\mathcal{L}(X)\). Alors, on s’intéresse au cas où

\[\mathrm{Var}[X - \alpha Y] \leq \mathrm{Var}[X]\]

Il est facile de vérifier qu’une condition nécessaire

\[\mathrm{Cov}[X,Y] \neq 0\]

Et si \(\alpha \in (0,\frac{-2\mathrm{Cov}[X,Y]}{\mathrm{Var}[Y]})\) (ou \(\alpha \in (\frac{-2\mathrm{Cov}[X,Y]}{\mathrm{Var}[Y]},0)\)), on a \(\mathrm{Var}[X - \alpha Y] \leq \mathrm{Var}[X]\).

De plus, \(\mathrm{Var}[X - \alpha Y]\) atteint son minimum lorsque \(\alpha = -\frac{\mathrm{Cov}(X,Y)}{\mathrm{Var}[Y]}\).

3.7 Algorithme de Metropolis-Hastings

Ici, on donne une version un peu plus générale que dans la feuille de TP : On utilise l’opérateur intégrale au lieu de la matrice de transitions. Le but est de comprendre pourquoi M-H nous donne un noyau de transition avec la bonne loi invariante et comment on fait la simulation. On remarque que tous les calculs seront identiques si on utilise la représentation matricielle. On va d’abord voir la construction de noyau de transition et vérifier qu’il nous donne la bonne loi invariante par les calculs algébriques. Ensuite, on donne deux décompositions pour mieux comprendre la dynamique associée et l’idée de la construction. On remarque que ce genre de dynamique nous permet de construire d’autres algorithmes pour des scénarios différents.

On suppose qu’on a une mesure \(\pi(dx)\) qui est absolument continue par rapport à la mesure (de Lebesgue) \(dx\), et un noyau de transition \(M(x,dy)\) tel que pour tous les \(x\), \(M(x,dy) = M(x,y)dy\) est une mesure de proba et est aussi absolument continue par rapport à \(dy\).

Alors, on définit le noyau Métropolis-Hastings par

\[ Q(x,dy) := r(x,y)M(x,dy) + \left(1 - \int r(x,s)M(x,ds) \right)\delta_x(dy) \] avec \[ r(x,y) = \frac{\pi(y)M(y,x)}{\pi(x)M(x,y)}\wedge 1 \]

Il est facile de vérifier que \[ \int \pi(dx)Q(x,dy)f(y) = \int f(y) \pi(dy) \] pour n’importe quelle fonction de test mesurable bornée \(f\).