Échantillonnage préférentiel

Introduction

Le problème est de calculer I(f) = \int f d\lambda,\mu est une mesure de probabilité sur un espace X et f : X \rightarrow \mathbb{R} est intégrable. Si \{X_n\} est une suite de variables aléatoires indépendantes et distribuées selon \mu, alors on peut approximer I(f) par

I_n(f) = \frac{1}{n}\sum_{i=1}^n f(X_i),

qui est dit un estimateur de Monte-Carlo.
En pratique, il peut être difficile de générer X_n \sim \lambda. On préférera alors introduire une mesure \mu, avec \lambda absolument continue par rapport à \mu, de sorte que

I_n(f;\mu) = \frac{1}{n}\sum_{i=1}^n f(Y_i) \tfrac{d\lambda}{d\mu}(Y_i), \quad Y_i \sim^{ind.} \mu,

soit une estimée de I(f) plus commode à calculer. Cette technique, dite de l’échantillonnage préférentiel, peut aussi servir à améliorer la qualité de l’estimateur I_n par exemple en réduisant sa variance.

Exemple.
Soit Z une variable aléatoire gaussienne de moyenne 0 et de variance 1, et considérons le problème d’estimer \mathbb{P}(Z > 10)= \int \mathbb{I}_{(10, \infty)} d\lambda\lambda est la distribution de Z. Une application directe de la méthode de Monte-Carlo serait d’utiliser I_n(f) = \frac{1}{n}\sum_{i=1}^n \mathbb{I}_{(\alpha, \infty)}(X_i), où X_i \sim^{ind.} \lambda. Le problème ici est qu’il faut prendre n de l’ordre de 10^{24} pour espérer que I_n(f) > 0. On peut autrement introduire \mu la distribution de |Z| + 10 pour obtenir

\mathbb{P}(Z > 10) = \int \frac{e^{50-10 x}}{2}e^{-(x+5)^2} d\mu(x).

Quelques trajectoires de I_n(f;\mu) sont tracées ci-dessous, avec P(Z > 10) notée en vert. Sur ces simulations, I_n(f) \equiv 0. Il faudrait prendre n de l’ordre de 10^{24} pour espérer que I_n(f) > 0.

echantillonnage-preferentiel-comparaison

 

Remarque.
L’intégration de Monte-Carlo et l’échantillonnage préférentiel sont particulièrement utile lorsque l’espace X, munit de la mesure de référence \mu, est trop complexe pour permettre l’application de méthodes numériques standard. Par exemple, X pourrait être l’espace des fonctions sur [0,1] et \mu la distribution d’un processus gaussien. On sait facilement échantillonner \mu pour produire \{x_i\} et I_n(f), mais il pourrait être difficile d’évaluer I(f) autrement.

Comparaison avec d’autres approches

On considère plus généralement l’approximation de l’intégrale I(f) par une moyenne

q_n(f) = \frac{1}{n} \sum_{i=1}^n f(x_i) \rho(x_i),

\{x_i\} est un ensemble de points et \rho une fonction de pondération.

Intégrale de Riemann-Stieltjes.
On utilise ici une grille déterministe de points \{x_i\}_{i=0}^n. Cette approche est surtout utilisée lorsque X = \mathbb{R} et \lambda est une mesure déterminée par la fonction de répartition F. Alors, pour une séquence de partitions P_n = \{x_{i,n}\}_{i=0}^n \subset [a,b] telle que \max_i F(x_{i,n})-F(x_{i-1,n}) \rightarrow 0, et si f est intégrable au sens de Riemann-Stieltjes, on a

\int_a^b f d\lambda = \lim_{n \rightarrow \infty} \sum_{i=1}^n f(x_{i,n})(F(x_{i,n}) - F(x_{i-1,n})).

Si F est continue, alors le choix x_{i,n} = F^{-1}(i/n) peut être considéré optimal: on trouve alors \int_a^b f d\lambda = \lim_n \frac{1}{n}\sum_{i=1}^n f(x_{i,n}). Si de plus \lambda est la mesure de Lebesgue et f continue, alors |q_n(f)-I(f)| = \mathcal{O}(n^{-1}). Par comparaison, en échantillonnage préférentiel, il est impossible d’obtenir une vitesse de convergence supérieure à o(n^{-1/2}), à moins que x\mapsto f(x)\rho(x) soit constante.

Analyse numérique bayésienne.
Le problème est ici probabilisé en quantifiant l’incertitude sur I(f) par le biais d’une distribution de probabilité sur f. On pourrait par exemple modéliser f par un processus gaussien ayant une moyenne m (disons m = 0) et une fonction de covariance C. La variance de l’estimateur q_n(f) = \frac{1}{n}\sum_{i=1}^n f(x_i) sera minimisée lorsque \sum_{i,j} C(x_i, x_j) est minimisée. Cette approche est particulièrement utile lorsqu’il est très couteux d’évaluer f en un point, et lorsqu’on a une idée approximative des endroits où fluctue le plus la fonction f. On prendra une fonction de covariance induisant moins de corrélation dans les zones de grande fluctuation de f et l’analyse numérique bayésienne propose une procédure systématique pour obtenir un choix optimal de points d’évaluation.

Propriétés des moyennes empiriques

Les estimées I_n(f) et I_n(f;\mu) s’expriment toutes deux sous la forme \overline X_n = \frac{1}{n}\sum_{i=1}^n X_i pour un certaine suite \{X_n\} de variables aléatoires i.i.d. telles que \mathbb{E}X_1 = I(f). La loi forte des grands nombres, sous l’hypothèse d’intégrabilité de X_n, nous assure que \overline X_n \rightarrow \mathbb{E}X_1, presque surement lorsque n \rightarrow \infty. On peut quantifier plus précisément la vitesse de cette convergence en termes de l’intégrabilité de |X_1|^p pour 2 > p \geq 1.

Théorème 1 (Marcinkiewicz-Zygmund).
Soit \{X_n\} une sequence de variables aléatoires i.i.d. telles que \mathbb{E}|X_1|^p < \infty pour un certain 1 \le p < 2. Alors,

|\overline X_n - {X_1}| = o\left(n^{1/p - 1}\right) \quad \text{p.s.}

Remarque.
On ne peut pas aller chercher une vitesse de convergence plus grande que 1/\sqrt{n}, même si disons \mathbb{E}{|X_n|^p} <\infty pour p \geq 2. C’est une conséquence du théorème central limite: dès que \text{Var}(X_n) < \infty, la distribution de \overline X_n est asymptotiquement déterminée.

Un comportement problématique des suites \{\overline X_k\}_{k\in \mathbb{N}} apparaissant dès que \mathbb{E}{|X_1|^p} = \infty, donne la réciproque au Théorème 1.

 

Proposition 2.
Soit \{X_n\} une séquence de variables aléatoires i.i.d. et intégrables Alors, pour tout M > 0, on a que

P\left( \left| \overline X_{n} - \overline X_{n-1}\right| > M n^{1/p-1} \;\; \text{ i.o.} \right)

est 0 ou 1 selon que \mathbb{E}{|X_1|^p} < \infty ou \mathbb{E}{|X_1|^p} = \infty.

Démonstration.
Sans perte de généralité, posons \mathbb{E}{X_n} = 0.

Supposons tout d’abord que \mathbb{E}{|X_n|^p} = \infty. Alors, avec probabilité 1, |X_n| > M n^{1/p} infiniment souvent pour tout M > 0. (En effet, \mathbb{E}{|X_n|^p} = \infty \Rightarrow \sum_n \mathbb{P}(|X_n| > Mn^{1/p}) = \infty et le second lemme de Borel-Cantelli nous donne \mathbb{P}(|X_n| > M n^{1/p} \text{ i.o.}) = 1 par l’indépendance des X_n.) De plus, comme X_n est intégrable, la loi forte des grands nombres implique \overline X_{n-1} \rightarrow 0 p.s. Comme |\overline X_n - \overline X_{n-1}| = |X_n - \overline X_{n-1}|/n \geq (|X_n| - |\overline X_{n-1}|)/n, on obtient le résultat annoncé.

Supposons maintenant que \mathbb{E}{|X_n|^p} < \infty. Alors, |X_n| \le M n^{1/p} p.s. éventuellement. Comme aussi |\overline X_{n-1}| \rightarrow 0 p.s., on en conclue que

|\overline X_n - \overline X_{n-1}| \le (|X_n| + |\overline X_{n-1}|)/n \le Mn^{1/p-1}

p.s. éventuellement. \Box

Corollaire 3.
Sous les hypothèses de la Proposition 2, si \mathbb{E}{|X_1|^p} = \infty alors |\overline X_n - \mathbb{E}{X_1}| \not = o\left(n^{1/p-1}\right) avec probabilité 1.

Bornes d’erreur locale

On s’intéresse ici à l’erreur relative

e_n(f;\mu) := \frac{|I(f) - I_n(f;\mu)|}{I(|f|)},

où on a normalisé par I(|f|) plutôt que I(f) par commodité mathématique.

Lemme 4 (voir Robert & Casella, 2004).
La variance de l’estimateur I_n(f; \mu) est minimisée en \mu = \lambda^*, où d \lambda^* / d\lambda = (|f|/\int|f| d\lambda)^{-1}.

Démonstration.
Notons que

\mathbb{E}{I_n(f;\mu)^2} \geq \frac{1}{n} I(|f|)^2

et

\mathbb{E}{I_n(f;\lambda^*)^2} = \frac{1}{n} \int \left(f \frac{d\lambda}{d\lambda^*}\right)^2 d \lambda^* = \frac{1}{n}I(|f|)^2.

Comme I_n(f;\mu) est sans biais, on en conclue que la borne inférieure \eqref{eq:borne_inf}, atteinte en \mu = \lambda^*, minimise la variance de I_n(f;\mu). $\Box$

On peut maintenant borner le carré de l’erreur relative en terme d’une divergence entre \lambda^* et \mu.

Proposition 5.
Soient \lambda^*, \mu, f et e_n(f; \mu) tels que précédemment. Alors,

\mathbb{E}{e_n(f;\mu)^2} \le \frac{D_{\chi^2}(\lambda^*, \mu) + 1}{n},

D_{\chi^2}(\lambda^*, \mu) = \int \left(\frac{d\mu}{d\lambda^*}-1\right)^2 d\mu est la divergence \chi^2 entre \lambda^* et \mu.

Démonstration.
On calcule simplement

\mathbb{E}{e_n(f;\mu)^2} = \frac{1}{n I(|f|)^2} \left( \int \left(f \frac{d\lambda}{d\mu} \right)^2 d\mu - I(f)^2\right)\\ \le \frac{1}{n} \int \left(\tfrac{|f|}{I(|f|)} \tfrac{d\lambda}{d\mu}\right)^2 d\mu\\ = \frac{1}{n} \int \left(\tfrac{d\lambda^*}{d\mu}\right)^2 d\mu = \frac{d_{\chi^2}(\lambda^*, \mu)+1}{n}.

\Box

Le prochain résultat permet de traiter le cas D_{\chi^2}(\lambda^*, \mu) = \infty \Leftrightarrow \text{Var}(I_n(f;\mu)) = \infty.

Proposition 6. (Chatterjee & Diaconis, 2015).
Soient \lambda^*, \mu, f et e_n(f; \mu) tels que précédemment. Alors, pour tout a > 0 et si X \sim \mu,

\mathbb{E}{e_n(f;\mu)} \le \sqrt{\frac{a}{n}} + 2 \sqrt{\mathbb{P}(\tfrac{d\mu}{d\lambda^*}(X) > a}).

Références

[1] Chatterjee, S. et Diaconis, P. (2017). The Sample Size Required in Importance Sampling. \textit{Arxiv preprint} (arXiv:1511.01437v3).

[2] Agapiou, S., Papaspiliopoulos, O., Sanz-Alonso, D. et Stuart, A. M. (2017). Importance Sampling: Intrinsic Dimension and Computational Cost. \textit{Arxiv preprint} (arXiv:1511.06196v3).

[3] Sanz-Alonso, D. (2016). Importance Sampling and Necessary Sample Size: an Information Theory Approach. \textit{Arxiv preprint} (arXiv:1608.08814v1).

[4] Robert, C. et Casella, G. (2004). \textit{Monte Carlo Statistical Methods}. Springer-Verlag.

[5] Diaconis, P. (1988). Bayesian Numerical Analysis. http://probabilistic-numerics.org/assets/pdf/Diaconis_1988.pdf

Leave a comment

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s