diff options
| -rw-r--r-- | Makefile | 2 | ||||
| -rw-r--r-- | doc/rapport.tex | 93 |
2 files changed, 49 insertions, 46 deletions
@@ -10,7 +10,7 @@ GSL_FLAGS:=$(shell pkg-config --libs gsl) NLOPT_FLAGS:=$(shell pkg-config --libs nlopt) VPATH = src:doc -.PHONY: clean +.PHONY: all clean all: rapport.pdf diff --git a/doc/rapport.tex b/doc/rapport.tex index 48ed0d7..41f1d84 100644 --- a/doc/rapport.tex +++ b/doc/rapport.tex @@ -13,14 +13,14 @@ \begin{document} \section*{Introduction} -Sujet : Comparer les performances des méthodes de stratification exposées dans \cite{etore:hal-00192540} +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 +On suppose qu'on cherche à calculer \begin{equation*} I=\int_{[0,1]^s}f(x)dx \end{equation*} @@ -30,20 +30,20 @@ par la formule : \begin{equation*} \lim_{N\rightarrow\infty}\frac{1}{N}\sum_{n=1}^Nf(\xi^{(n)}) = 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 +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*} +\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. +\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*} @@ -53,11 +53,11 @@ 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 : +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 : +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} @@ -71,11 +71,12 @@ asymptotique usuel pour $\mu$ de la forme : $[\bar{X}_N-c_{\alpha}\frac{\sigma}{ avec $\alpha$ le niveau désiré pour l'intervalle et $c_{\alpha}$ le quantile en $1-\frac{\alpha}{2}$ d'une loi gaussienne centrée réduite. Dans l'article \cite{tuffin2004randomization}, Tuffin suggère d'utiliser la borne de Berry-Esseen pour le cas non asymptotique. Si on introduit les quantités -$\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 +$\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}} +\forall x, \quad |F_n(x) - \mathcal{N}(x)| \leq C\frac{\beta}{\sigma^3 \sqrt{N}} \end{equation*} +pour une constante universelle $C$\footnote{La meilleure valeur connue à l'heure actuelle est $C< 0.4748$.}. 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*} @@ -99,14 +100,14 @@ Alors que Tuffin trouve : ce qui me semble être une erreur. Si on veut comparer les demi-tailles des intervalles de confiance, on voit que l'intervalle de Tuffin rajoute $\frac{C\beta}{\sigma^2N}$ à l'intervalle habituel assymptotique tandis que je trouve un intervalle possiblement égal à $\mathbb{R}$ tout entier pour des petites valeurs de $N$. -C'est pourquoi dans les applications ultérieures j'utilise l'intervalle classique $[\bar{X}_N-c_{\alpha}\frac{\sigma}{\sqrt{N}}, +C'est pourquoi dans les applications ultérieures, j'utilise l'intervalle classique $[\bar{X}_N-c_{\alpha}\frac{\sigma}{\sqrt{N}}, \bar{X}_N+c_{\alpha}\frac{\alpha}{\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 +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}))$. @@ -116,10 +117,10 @@ $\frac{1}{N_{i}}\sum_{j=1}^{N_{i}}f(X_{i}^{j})$ où $N_{i}$ est le nombre de tir \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. +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. +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}}= @@ -127,83 +128,84 @@ Si on note $\sigma_{i} = \mathrm{Var}(f(X_{i})) = \mathrm{Var}(f(X)|X\in A_{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}) = +\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( +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 +\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*} + +\begin{multline} +\label{eq:borninf} +\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=\frac{1}{N}(\sum_{i=1}^Ip_i\sigma_i)^2=\frac{1}{N}\sigma_\star^2 +\end{multline} + Et si \begin{equation} \label{eqqi} -\forall i \in \{1,\ldots,I\} \quad q_i=\frac{p_i\sigma_i}{\sum_{j=1}^Ip_j\sigma_j} +\forall i \in \{1,\ldots,I\}, \quad q_i=\frac{p_i\sigma_i}{\sum_{j=1}^Ip_j\sigma_j} \end{equation} -$\textrm{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é, +$\textrm{Var}(\mu_{strat})$ atteindra la borne inférieure $\sigma_\star^2/N$. 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 +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 +en choisissant une répartition déterministe, par exemple l'allocation proportionelle $q_i=p_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. -Essayons maintenant de décrire plus précisément comment calculer précisément la proportion de tirages à allouer par strates. En effet, il y a deux contraintes +Essayons maintenant de décrire plus précisément comment calculer la proportion de tirages à allouer par strates. En effet, en pratique il y a deux contraintes en plus de la formule qui nous donne la suite des $q_i$. Tout d'abord le nombre de tirages doit évidemment être entier, et de plus il doit y avoir au moins -un tirage par strate, car Étoré et Jourdain prouvent qu'il y a convergence de leurs estimateurs de variances conditionnelles $\hat{\sigma}_i$ si il y a au moins un -tirages effectué par strates à chaque étape. -Notons $I$ le nombre de strates, $\hat{\sigma}_i^k$ l'estimateur à l'étape $k$ de la variance de $X$ sachant que $X$ est dans le $i_{\textrm{ème}}$ intervalle -, $N_i^k$ le nombre total de tirages effectués de l'étape 0 à l'étape $k$ dans la $i_{\textrm{ème}}$ strate et $M_i^k$ le nombre de tirages à effectuer à +un tirage par strate, sinon on ne peut garantir la convergence des variances conditionelles empiriques $\hat{\sigma}_i$ vers la vraie valeur. +Notons $I$ le nombre de strates, $\hat{\sigma}_i^k$ l'estimateur à l'étape $k$ de la variance de $X$ sachant que $X \in A_i$, $N_i^k$ le nombre total de tirages effectués de l'étape 0 à l'étape $k$ dans la $i^{\textrm{ème}}$ strate et $M_i^k$ le nombre de tirages à effectuer à l'étape $k$ dans la strate $i$ qui est donc $M_i^k = N_i^k - N_i^{k-1}$. -Étape 1 : +\paragraph{Étape 1 :} \begin{equation*} -\forall i \in I \quad M_i^1 = N_1/I +\forall i \in I, \quad M_i^1 = 1+p_iN_1 \end{equation*} - -Étape k : -On commence d'abord par calculer les $\hat{\sigma}_i^{k-1}$ avec les tirages de l'étape précédente. +C'est l'allocation proportionelle plus un, pour garantir au moins un tirage par strate. Si les $M_i^1$ ne sont pas entiers, voir la méthode ci-dessous pour les arrondir. +\paragraph{Étape k :\\} +On commence d'abord par calculer les $\hat{\sigma}_i^{k-1}$ avec les tirages de l'étape précédente. \begin{equation*} \forall i \in I \quad \hat{\sigma}_i^{k-1} = \sqrt{\frac{1}{N_i^{k-1}}(\sum_{j=1}^{N_i^{k-1}}(f(X_i^j))^2-(\frac{1}{N_i^{k-1}}\sum_{j=1}^{N_i^{k-1}}(f(X_i^j))^2)} \end{equation*} Ensuite on peut calculer les $M_i^k$. On sait qu'il nous faut au moins un tirage par strate on peut donc écrire $M_i^k$ de sous la forme $m_i^k + 1$ avec -$m_i^k \in \mathbb{N}$ et comme on a $M_i^k = N_i^k - N_i^{k-1}$ il vient $\sum_{i=1}^Im_i^k=N^k-N^{k-1}-I$. Les auteurs de l'article proposent deux méthodes +$m_i^k \in \mathbb{N}$ et comme on a $M_i^k = N_i^k - N_i^{k-1}$ il vient $\sum_{i=1}^Im_i^k=N^k-N^{k-1}-I$. Les auteurs de l'article proposent deux méthodes pour calculer les $m_i^k$, je n'en ai utlisé qu'une seule que je vais présenter ici. -D'après \eqref{eqqi} on connait la proportion optimale à faire par strates en fonction des $\hat{\sigma}_i^{k-1}$. Sachant de plus que $\sum_{i=1}^Iq_i=1$, +D'après \eqref{eqqi}, on connait la proportion optimale à faire par strates en fonction des $\hat{\sigma}_i^{k-1}$. Sachant de plus que $\sum_{i=1}^Iq_i=1$, on voit qu'un choix de $m_i^k$ idéal serait : \begin{equation*} -\forall i \in I \quad (m_i^k)^* = \frac{p_i\hat{\sigma}_i^{k-1}}{\sum_{j=1}^Ip_j\hat{\sigma}_j^{k-1}}(N^k-N^{k-1}-I) +\forall i \in I, \quad (m_i^k)^* = \frac{p_i\hat{\sigma}_i^{k-1}}{\sum_{j=1}^Ip_j\hat{\sigma}_j^{k-1}}(N^k-N^{k-1}-I) \end{equation*} -Pour que les $m_i^k$ soit les entiers au plus près de leurs valeurs optimales tout en assurant que leur somme ne dépasse pas le nombre total de tirages à faire -à l'étape k. On pose $m_1^k =\floor{(m_1^k)^*}$. Puis pour chaque $i>1$ on pose : +Pour que les $m_i^k$ soient les entiers au plus près de leurs valeurs optimales, tout en assurant que leur somme ne dépasse pas le nombre total de tirages à faire +à l'étape $k$, Étoré et Jourdain proposent la méthode suivante. On pose d'abort $m_1^k =\floor{(m_1^k)^*}$. Puis pour chaque $i>1$ on pose : \begin{equation*} m_i^k=\floor{(m_1^k)^*+\ldots+(m_i^k)^*}-\floor{(m_1^k)^*+\ldots+(m_{i-1}^k)^*} \end{equation*} Ce choix des $m_i^k$ assure que notre nombre de tirages par strates sera à une distance d'au plus $1$ de l'allocation optimale. Étoré et Jourdain prouvent enfin le résultat suivant qui est un théorème centrale limite pour la loi limite de leur estimateur $\mu_{strat}$. -Si $\E[f^2(X)]<\infty$ et si $k/N ^k\to0$ quand $k\to0$ alors : +Si $\E[f^2(X)]<\infty$ et si $k/N ^k\to0$ quand $k\to \infty$ alors : \begin{equation*} -\sqrt{N^k}(\mu_{strat}^k-\mu)\overset{en loi}{\to}\mathcal{N}(0,\sigma^2_*) +\sqrt{N^k}(\mu_{strat}^k-\mu)\underset{k\to\infty}{\overset{d}{\longrightarrow}}\mathcal{N}(0,\sigma^2_*) \end{equation*} -où $\sigma^2_*$ désigne la borne inférieur de $\textrm{Var}(\mu_{strat})$ trouvée plus haut. +avec $\sigma_*=\sum_i^Ip_i\sigma_i$, ce qui atteint la borne inférieure de $\textrm{Var}(\mu_{strat})$ trouvée en \eqref{eq:borninf}. \section{Exemples} Nous allons maintenant implémenter et comparer les deux techniques d'estimateur d'espérance sur deux exemples décrit dans \cite{etore:hal-00192540} \subsection{Calcul de l'espérance d'une loi normale} -Dans ce premier exemple assez simple, on va chercher à calculer $ \mu = \E[X]$ où $X$ suit une loi normale centrée réduite +Dans ce premier exemple assez simple, on va chercher à calculer $ \mu = \E[X]$ où $X$ suit une loi normale centrée réduite \subsubsection{Échantillonnage stratifié} Pour cet exemple, les strates $A_i$ seront au nombre de $10$ et \[\forall i \in{\{1,\ldots,10\}}\quad A_i =]y_{i-1}, y_i] @@ -304,6 +306,7 @@ Pour conculure on peut constater que quasi Monte-Carlo randomisé semble être t cas sur les deux exemples étudiés ici. Sachant de plus que quasi Monte-Carlo randomisé est bien plus simple à implémenter, l'intérêt pratique de l'échantillonage stratifié pour étudier des cas simples ne semble pas démontré. +\newpage \printbibliography \end{document} |