\documentclass[a4paper, 11pt, french]{article}
\usepackage{babel}
\usepackage{fontspec}
\usepackage{amssymb, amsmath}
\usepackage[pagebackref=true,breaklinks=true,colorlinks=true,backref=false]{hyperref}
\usepackage[capitalize,noabbrev]{cleveref}
\usepackage[backend=biber,style=trad-plain]{biblatex}
\addbibresource{rapport.bib}
\newcommand{\E}{\mathbb{E}}
\newcommand{\V}{\mathbb{V}}

\begin{document}
\section*{Introduction}
Sujet : Comparer les performances des méthodes de stratification exposées dans \cite{etore:hal-00192540} 
et de QMC randomisées exposées dans \cite{tuffin2004randomization} sur les exemples présentés dans \cite{etore:hal-00192540}

Références : \cite{etore:hal-00192540} \cite{tuffin2004randomization}

\section{Quasi Monte Carlo randomisé}
\subsection{Présentation mathématique}
On suppose qu'on cherche à calculer  
\begin{equation*}
I=\int_{[0,1]^s}f(x)dx
\end{equation*}
Soit $(\xi^{(n)})_{n\in\mathbb{N}}$ une suite à discrépance faible à valeurs dans $[0,1]^s$.
Au bout de N tirages de cette suite, la méthode de quasi Monte Carlo  nous donnerait une approximation de la valeur de $I$
par la formule :
\begin{equation*}
\frac{1}{N}\sum_{n=1}^Nf(\xi^{(n)})\overset{p.s}{\to} I
\end{equation*}
Le souci de cette méthode est de ne pas pouvoir obtenir d'erreur de l'estimateur facilement. En effet l'inégalité de Koksma–Hlawka nous donne une borne de 
cette erreur mais nécessite de pouvoir calculer la variation finie de la fonction $f$ ce qui n'est pas toujours possible en pratique. C'est pourquoi on introduit
la méthode de quasi Monte Carlo randomisé. Elle consiste à rajouter un décalage aléatoire à la suite $(\xi^{(n)})$. Soit $X$ une variable aléatoire
uniformément distribuée sur ${[0,1]^s}$. On s'intéresse dorénavant à la somme aléatoire :
\begin{equation*}
Z=\frac{1}{N}\sum_{n=1}^Nf(\{\xi^{(n)}+X\})
\end{equation*} 
où \{$X$\} désigne la partie fractionnaire de $X$.
L'idée consiste à faire $K$ tirages de $Z$, ce qui revient à faire $K$ tirages de $X$ vu que la suite $(\xi^{(n)})$ est déterministe, puis de faire
une méthode de Monte Carlo classique sur ces $K$ tirages. On a donc :
\begin{equation}
\label{eqrqmc}\frac{1}{K}\sum_{k=1}^KZ_{k}=\frac{1}{K}\sum_{k=1}^K\frac{1}{N}\sum_{n=1}^Nf(\{\xi^{(n)}+X_{k}\})\overset{p.s}{\to} I
\end{equation} 
Il convient maintenant d'estimer l'avantage de cette méthode par rapport à une technique de Monte-Carlo usuelle, en terme de réduction de variance. 
On voit dans \eqref{eqrqmc} qu'il nous coûte $NK$ évaluations de la fonction $f$ pour calculer notre estimation de $I$.
Donc la variance pour la méthode de Monte-Carlo de base serait :
\begin{equation*}
\mathrm{Var}(\frac{1}{NK}\sum_{n=1}^{NK}f(X_n))=\frac{1}{NK}\mathrm{Var}(f(X))
\end{equation*}
Il nous faut la comparer avec :
\begin{equation*}
\mathrm{Var}(\frac{1}{K}\sum_{k=1}^K\frac{1}{N}\sum_{n=1}^Nf(\{\xi^{(n)}+X_{k}\}))=\frac{1}{K}\mathrm{Var}(\frac{1}{N}\sum_{n=1}^Nf(\{\xi^{(n)}+X\}))
\end{equation*}
On a donc réduction de variance si et seulement si : 
\begin{equation*}
\mathrm{Var}(\frac{1}{N}\sum_{n=1}^Nf(\{\xi^{(n)}+X\}))<\frac{1}{N}\mathrm{Var}(f(X))
\end{equation*}
Or dans \cite{Tuffin:1997:VRA:268403.268419} il est montré que pour $f$ une fonction à variations bornées on a : 
\begin{equation}
\mathrm{Var}(\frac{1}{N}\sum_{n=1}^Nf(\{\xi^{(n)}+X\}))=O(N^{-2}(\log N)^{2s})
\end{equation}
Pour cette classe de fonctions on a une réduction de variance par un facteur de $N / (\log N)^{2s}$.

\subsection{Calcul de vitesse et d'intervalle de confiance}
On s'interesse maintenant à la convergence vers la loi normale de notre estimateur. On rappelle d'abord le cas d'une méthode de Monte-Carlo classique. 
Soit $X$ une variable aléatoire, on note $\mu = \E[f(X)]$ son espérance qu'on cherche à calculer et $\sigma$ son écart type. Alors si les $(X_n)_{1\leq n\leq N}$ 
sont $N$ tirages de la variable $X$, $\bar{X}_N=\frac{1}{N}\sum_{n=1}^Nf(X_n)$ sera notre estimateur de $\mu$ et on aura l'intervalle de confiance 
assymptotique usuel pour $\mu$ de la forme : $\mathopen{[}\bar{X}_N-c_{\alpha}\frac{\sigma}{\sqrt{N}};\bar{X}_N+c_{\sigma}\frac{\alpha}{\sqrt{N}}\mathclose{]}$ 
avec $\alpha$ le niveau désiré pour l'intervalle et $c_{\alpha}$ le quantile de $1-\frac{\alpha}{2}$ d'une loi gaussienne centrée réduite.

Dans son article \cite{tuffin2004randomization} M.Tuffin suggère d'utiliser la borne de Berry-Esseen pour le cas non assymptotique. Si on introduit la variable 
$\beta = \E[|f(X)-\E[f(X)]|^3]$, $F_n$ la fonction de répartition de $\sqrt{N}\frac{\bar{X}_N-\mu}{\sigma}$et $\mathcal{N}$ la fonction de répartition d'une
loi normale centrée réduite, alors l'inégalité de Berry-Esseen affirme que :
\begin{equation*}
\forall x \quad |F_n(x) - \mathcal{N}(x)| \leq C\frac{\beta}{\sigma ^3 \sqrt{N}}
\end{equation*}
Ce qui nous donne $\mathcal{N}(x) - C\frac{\beta}{\sigma ^3 \sqrt{N}} \leq F_n(x) \leq \mathcal{N}(x) + C\frac{\beta}{\sigma ^3 \sqrt{N}} $ et en 
soustrayant l'encadrement équivalent de $F_n(-x)$ il vient :
\begin{equation*}
F_n(x)-F_n(-x) \geq \mathcal{N}(x) - \mathcal{N}(-x) -2C\frac{\beta}{\sigma ^3 \sqrt{N}}
\end{equation*}
On veut maintenant égaliser cette borne inférieure avec $1- \alpha$ pour trouver les bornes de l'intervalle de confiance.
Sachant que $\mathcal{N}(x) = 1 - \mathcal{N}(-x)$ on a $1 - \alpha =  1 - 2\mathcal{N}(-x) - 2C\frac{\beta}{\sigma ^3 \sqrt{N}}$ et on trouve en fin de 
compte :
 $x = -\mathcal{N}^{-1}(\frac{\alpha}{2} - C\frac{\beta}{\sigma ^3 \sqrt{N}})$.

\section{Échantillonnage stratifié}
\subsection{Présentation mathématique}
Soit $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ et $X$ un vecteur aléatoire à $d$ dimensions.
On cherche à calculer $\mu=\E(f(X))$.
On suppose maintenant qu'on divise $\mathbb{R}^{d}$ en $I$ strates appelées $(A_{i})_{1\le i \le I}$ et de telle sorte que 
$\forall i \in {1,...,i}$ on ait $\mathbb{P}(X\in A_{i}) > 0$ . On notera ces probabilités $p_{i}$.
Si on note maintenant $(X_{i})_{1\le i \le I}$ l'ensemble de variables aléatoires suivant la loi conditionnelle de $X$ sachant $X\in A_{i}$,
on a clairement $\E(f(X)) = \sum_{i=1}^{I}p_{i}\E(f(X_{i}))$.
Or si on sait simuler chacune des variables $(X_{i})$, on peut estimer chacune des espérances $\E(f(X_{i}))$ par un estimateur de Monte-Carlo usuel
$\frac{1}{N_{i}}\sum_{j=1}^{N_{i}}f(X_{i}^{j})$ où $N_{i}$ est le nombre de tirages de la variable $X_{i}$. Et donc il vient naturellement :
\begin{equation*}
\mu_{strat}=\sum_{i=1}^{I}p_{i}\sum_{j=1}^{N_{i}}f(X_{i}^{j})=
\frac{1}{N}\sum_{i=1}^{I}\frac{p_{i}}{q_{i}}\sum_{j=1}^{q_{i}N}f(X_{i}^{j}) \overset{p.s}{\to} \mu
\end{equation*}
où $N$ désigne le nombre total de tirages, c'est à dire $N=\sum_{i=1}^{I}N_{i}$, et $q_{i}=N_i / N$ la proportion de tirages dans chacune 
des strates. 

On veut maintenant regarder dans quels cas cette méthode peut être utile pour réduire la variance. 
Si on note $\sigma_{i} = \mathrm{Var}(f(X_{i})) = \mathrm{Var}(f(X)|X\in A_{i})$, la variance de notre estimateur est :
\begin{equation*}
\mathrm{Var}(\mu_{strat})=\sum_{i=1}^{I}\frac{p_{i}^{2}\sigma_{i}^2}{N_{i}}=
\frac{1}{N}\sum_{i=1}^{I}\frac{p_{i}^{2}\sigma_{i}^2}{q_{i}}
\end{equation*}
Tandis qu'avec un estimateur usuel de Monte-Carlo $\frac{1}{N}\sum_{j=1}^{N}f(X^{j})$, noté $\mu_{mc}$, pour le même nombre $N$ total de tirages, on aurait :
\begin{equation*}
\mathrm{Var}(\mu_{mc}) = 
\frac{1}{N}\left(\sum_{i=1}^{I}p_{i}(\sigma_{i}^2 + \E^{2}[f(X_{i})])-(\sum_{i=1}^{I}p_{i}\E[f(X_{i})])^{2}\right)
 \geq \frac{1}{N}\sum_{i=1}^{I}p_i\sigma_i^2
\end{equation*}
où la dernière inégalité est due au fait que par l'inégalité de Jensen appliquée à $\E[f(X_{i})]$ avec la suite des $p_i$ comme probabilité discrète(
$\sum_{i=1}^{I}p_i = 1$) on a :
\begin{equation*}
\frac{1}{N}\sum_{i=1}^{I}p_i\E^{2}[f(X_{i})] \geq (\frac{1}{N}\sum_{i=1}^{I}p_i\E[f(X_{i})])^2
\end{equation*}
Si on pose $q_i = p_i$, $\mathrm{Var}(\mu_{strat})$ atteindra la borne inférieure trouvée pour $\mathrm{Var}(\mu_{mc})$.
Mais on peut trouver un choix de $q_i$ plus optimal en terme de réduction de variance.
En effet, toujours par l'inégalité de Jensen avec cette fois la suite des $q_i$ vue comme une probabilité discrète, on a :
\begin{equation*}
\mathrm{Var}(\mu_{strat})=\frac{1}{N}\sum_{i=1}^{I}\frac{p_i^2\sigma_i^2}{q_i}=\frac{1}{N}\sum_{i=1}^I(\frac{p_i\sigma_i}{q_i})^2q_i
\geq\frac{1}{N}(\sum_{i=1}^I\frac{p_i\sigma_i}{q_i}q_i)^2
\end{equation*}
Et si
\begin{equation}
\label{eqqi}
\forall i \in \{1,...,I\} \quad q_i=\frac{p_i\sigma_i}{\sum_{j=1}^Ip_j\sigma_j}
\end{equation}
${Var}(\mu_{strat})$ atteindra la borne inférieure trouvée précédemment. Maintenant qu'on a décrit le principe de la méthode de l'échantillonnage stratifié,
et qu'on a trouvé quelle est la proportion optimale de tirages à faire par strates, intéressons nous à l'algorithme pour implémenter ce modèle.

\subsection{Description de l'algorithme}
On a vu dans la section précédente que les $q_i$ optimaux dépendent de la suite $(\sigma_i)_{1\leq i\leq I}$ qui sont les variances conditionnelles de $X$
sachant $X\in A_{i}$. Mais si ces variances sont inconnues il nous faut les estimer. On ne peut trouver la répartition optimale de tirages à faire par strates 
du premier coup. On utilisera donc un algorithme adaptatif. On commence par tirer un nombre $N_1$ assez restreint de variables $(X_{i})_{1\leq i \leq I}$
en choisissant une répartition déterministe, par exemple la répartition uniforme $q_i=\frac{N_1}{I}$ ce qui nous donnera un premier estimateur de  $\mu=\E(f(X))$ mais
surtout les premiers estimateurs de $\sigma_{i}$. Ensuite on pourra tirer un nombre $N_2 > N_1$ de $(X_{i})$ avec de nouveaux $q_i$ donnés par \eqref{eqqi}
ce qui nous donnera un nouvel estimateur de $\mu$ (en utilisant les $N_1 + N_2$ premières valeurs) puis de nouveaux $q_i$ plus précis, et ainsi de suite.


\input{table.tex}

\printbibliography

\end{document}