{ "cells": [ { "cell_type": "markdown", "id": "7a149ebe-0121-407d-ab82-22dbdae6833e", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "# Réduction de variance: préconditionnement et importance sampling\n", "\n", "Sur on problème de modélisation classique en assurance, on illustre l'importance de l'erreur relative puis deux techniques de réduction de variance: \n", "\n", "- méthode par préconditionnement, \n", "- échantillonage d'importance (important pour les événements rares)." ] }, { "cell_type": "code", "execution_count": 1, "id": "1bb383e8-5124-4169-b672-77bd10f1b6d1", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy import stats\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "sns.set_theme() \n", "from numpy.random import default_rng\n", "rng = default_rng()" ] }, { "cell_type": "markdown", "id": "e93be7c7-b92f-44c1-9dde-c4e06aaa199c", "metadata": {}, "source": [ "## Charge sinistre et loi Poisson-composée\n", "\n", "On définit la _charge sinistre totale_ (sur une période $T$) par la variable aléatoire positive\n", "\\begin{equation*}\n", " S = \\sum_{i = 1}^N X_i\n", "\\end{equation*}\n", "où $N$ est une variable aléatoire à valeurs dans $\\mathbf{N}$ représentant le nombre de sinistres sur la période $T$, et pour $i \\ge 1$, $X_i$ est une variable aléatoire à valeurs dans $\\mathbf{R}_+$ représentant le coût du i-ème sinistre, avec la convention selon laquelle la somme est nulle si $N = 0$. Les $(X_i)_{i \\ge 1}$ sont supposées indépendantes et identiquement distribuées, et indépendantes de $N$ (indépendance fréquences - coûts).\n", "\n", "Une modélisation classique est de considérer \n", "\n", "- $N$ de loi de Poisson de paramètre $\\lambda > 0$, \n", "- $X_1$ de loi log-normale de paramètres $\\mu > 0$, $\\sigma^2 > 0$, c'est à dire $X_1 = \\exp(G_1)$ avec $G_1 \\sim \\mathcal{N}(\\mu, \\sigma^2)$. \n", "\n", "Le but est d'estimer la **probabilité de dépassement** c'est à dire calculer la probabilité que la charge sinistre totale dépasse un seuil $K$:\n", "\\begin{equation*}\n", " p = \\mathbf{P}\\bigl[S > K\\bigr] \\quad \\text{pour $K$ grand}\n", "\\end{equation*}\n", "\n", "Dans la suite on prend $\\lambda = 10$, $\\mu = 0.1$ et $\\sigma = 0.3$ et on considère plusieurs valeurs du seuil $K$." ] }, { "cell_type": "markdown", "id": "6851b921-b5ae-45b0-8f2c-23b9590a0b7b", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "question" ] }, "source": [ "### Question: simulation de la charge sinistre totale \n", "\n", "Ecrire une fonction `simu_S(size, mu, sigma, lambd)` qui renvoie un échantillon de taille `size` de réalisations indépendantes de $S$." ] }, { "cell_type": "code", "execution_count": null, "id": "c7d80c68-990f-4101-8a52-44e0fe73d476", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "aremplir" ] }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "53745a89-b51c-480c-ae1c-f2a5a3f09f0f", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "question" ] }, "source": [ "### Question: représentation graphique \n", "\n", "Représenter l'histogramme d'un échantillon de $100\\,000$ réalisations de $S$ et du seuil $K = 20$ par une ligne verticale rouge." ] }, { "cell_type": "code", "execution_count": null, "id": "6842db6d-f50d-4c77-bc5f-6caf48ac8354", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "aremplir" ] }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "988dc562-8060-4a08-afb2-98396ee4a032", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## Estimateur Monte Carlo et erreur relative \n", "\n", "Soit $p_n = \\displaystyle \\frac{1}{n} \\sum_{j=1}^n \\mathbf{1}_{S^{(j)} > K}$ l'estimateur Monte Carlo de $p= \\mathbf{P}\\bigl[S > K\\bigr]$ où $(S^{(j)})_{j=1,\\dots,n}$ est une suite _i.i.d._ de même loi que $S$.\n", "\n", "On rappelle que:\n", "\n", "- l'**erreur absolue** de l'estimateur Monte Carlo $p_n$ est définie par $|p_n - p|$ et qu'avec probabilité 0.95 cette erreur est bornée par $e_n = 1.96 \\frac{\\sigma_n}{\\sqrt{n}}$ avec $\\sigma_n^2 = p_n - p_n^2$,\n", "- l'**erreur relative** de l'estimateur Monte Carlo est définie par $\\frac{|p_n - p|}{p}$ que l'on majore avec probabilité 0.95 par $\\frac{e_n}{p_n}$." ] }, { "cell_type": "markdown", "id": "50ea185f-b1b6-450d-84ce-8ea6c2964654", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "question" ] }, "source": [ "### Question: erreur relative \n", "\n", "Ecrire une fonction `relative_error` qui à partir d'un échantillon de $S$ (de taille $n$) et d'une valeur de seuil $K$ renvoie la probabilité $p_n$ et l'erreur relative (plus exactement la borne $\\frac{e_n}{p_n}$ à 95%). \n", "\n", "Tracer l'erreur relative d'un échantillon de taille $100\\,000$ en fonction de $K$ pour $K$ allant de 20 à 30. Comment interpréter cette courbe? " ] }, { "cell_type": "code", "execution_count": null, "id": "3e487ddc-f33f-400b-a03f-dfe7388e1afa", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "aremplir" ] }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "ee4c5ab5-64f7-4495-a97b-7da7c60fcfc3", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "question" ] }, "source": [ "### Question: Monte Carlo à précision fixée\n", "\n", "Mettre en oeuvre un estimateur de Monte Carlo qui s'arrête dès que l'erreur relative est de 5%. On pourra par exemple introduire la variable aléatoire\n", "\\begin{equation*}\n", " \\tau^{(m)} = \\inf \\bigl\\{n \\ge 1, e_{n m} \\le 0.05 p_{n m} \\bigr\\},\n", "\\end{equation*}\n", "qui dépend d'un paramètre $m$ fixé, par exemple $m = 10\\,000$, et renvoyer $p_{\\tau^{(m)}}$ ainsi que l'erreur relative et la taille de l'estimateur associé. \n", "Le paramètre $m$ permet de recalculer l'estimateur et l'erreur uniquement toutes les $m$ itérations et donc de réduire la complexité par rapport au choix naïf $m = 1$. On appelle ce paramètre $m$ la taille du _batch_ (`batch_size`). \n", "Le nombre d'itérations (la taille de l'échantillon) dans la méthode de Monte Carlo pour un $\\tau^{(m)}$ donné est donc $\\tau^{(m)} \\times m$. \n", "\n", "Définir la fonction qui code cet estimateur Monte Carlo:\n", "```\n", "monte_carlo_relative(mu, sigma, lambd, K, batch_size = 10000, error = 0.05)\n", "```" ] }, { "cell_type": "code", "execution_count": null, "id": "74ff1d4b-d9bf-4e1a-977f-3cd862e698a9", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "aremplir" ] }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "5dde5182-f412-4166-93a3-c896cdd0db9f", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "question" ] }, "source": [ "### Question: complexité en fonction de $K$\n", "\n", "Reproduire un tableau de résultat similaire au tableau suivant obtenu avec cet estimateur de Monte Carlo adaptatif jusqu'à l'itération $\\tau^{(m)} \\times m$ pour une erreur relative de 10% et pour différentes valeurs de $K = 20,\\dots,25$. \n", "\n", "Tracer le nombre d'itérations nécessaires en fonction de $K$." ] }, { "cell_type": "code", "execution_count": 12, "id": "2067ee08-df59-44e1-9699-38f57d8cd73a", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "question" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Probabilité $p_n$Erreur relativeItérations
200.0216000.09327720000
210.0131000.09821930000
220.0078200.09873350000
230.0039100.098927100000
240.0022590.099908170000
250.0012580.099186310000
\n", "
" ], "text/plain": [ " Probabilité $p_n$ Erreur relative Itérations\n", "20 0.021600 0.093277 20000\n", "21 0.013100 0.098219 30000\n", "22 0.007820 0.098733 50000\n", "23 0.003910 0.098927 100000\n", "24 0.002259 0.099908 170000\n", "25 0.001258 0.099186 310000" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import pandas as pd\n", "df = pd.read_pickle(\"data/iterations_df.pkl\")\n", "df" ] }, { "cell_type": "code", "execution_count": null, "id": "92e26380-e252-4d58-bb52-5c0e28660cbd", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "aremplir" ] }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "9eb309d4-44f9-430c-9cdd-6c58489abf9c", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## Réduction de variance par préconditionnement \n", "\n", "Pour réduire la variance on utilise d'abord une méthode de préconditionnement, c'est à dire qu'on considère la variable aléatoire \n", "\\begin{equation*}\n", " M = \\inf \\bigl\\{ r \\ge 1, \\sum_{i=1}^r X_i > K \\bigr\\}\n", "\\end{equation*}\n", "et la représentation suivante\n", "\\begin{equation*}\n", " p = \\mathbf{E} \\bigl[\\phi(M) \\bigr] \\quad \\text{avec} \\quad \\phi(m) = \\mathbf{P} \\bigl[ N \\ge m \\bigr].\n", "\\end{equation*}" ] }, { "cell_type": "markdown", "id": "8ef32c82-ad82-4087-ac8d-4bb60c325cf3", "metadata": { "tags": [ "question" ] }, "source": [ "### Question: simulation de $M$ \n", "\n", "Ecrire une fonction `simu_M` similaire à la fonction `simu_S` avec l'argument $K$ supplémentaire qui renvoie un échantillon _i.i.d._ de même loi que $M$. " ] }, { "cell_type": "code", "execution_count": null, "id": "016580ba-e210-46ef-925b-bcb7f4221dc1", "metadata": { "tags": [ "aremplir" ] }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "e7d2005e-4b92-4796-a47f-0eeda856f72f", "metadata": { "tags": [ "question" ] }, "source": [ "### Question: Monte Carlo et ratio de variance \n", "\n", "En utilisant la fonction `monte_carlo` du TP précédent. Calculer le ratio de variance entre l'estimateur $p_n$ et l'estimateur basé sur la représentation $p = \\mathbf{E}[\\phi(M)]$ où $\\phi$ est calculée en utilisant la fonction de survie et la fonction de masse de la loi de Poisson (cf. la documentation de `stats.poisson`). Faire ce calcul pour différentes valeurs de $K$ et $n = 20\\,000$" ] }, { "cell_type": "code", "execution_count": null, "id": "986219d3-c538-4c1e-87b5-56d1c80d15aa", "metadata": { "tags": [ "aremplir" ] }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "9d2ad516-a1c2-4b07-b0a9-bb254ba5f684", "metadata": {}, "source": [ "## Réduction de variance par échantillonage d'importance\n", "\n", "Pour réduire la variance sans faire exploser la complexité pour les grandes valeurs de $K$ on propose une méthode d'échantillonage d'importance (Importance Sampling) en modifiant la loi de la variable aléatoire $N$ (on peut faire un autre choix, en changeant la loi des $X_i$ ou bien en changeant la loi de $N$ et des $X_i$). Le changement de loi proposé ici repose sur le changement de probabilité, pour $\\theta \\in \\mathbf{R}$\n", "\\begin{equation*}\n", " \\frac{\\mathrm{d} \\mathbf{P}}{\\mathrm{d} \\mathbf{P}_\\theta} = L_\\theta \\quad \\text{avec} \\quad L_\\theta = \\exp\\bigl(-\\theta N + \\psi(\\theta) \\bigr),\n", "\\end{equation*}\n", "où $\\psi(\\theta) = \\log \\mathbf{E} \\bigl[ \\exp(\\theta N) \\bigr] = \\lambda (e^\\theta - 1)$. On vérifie par le calcul que la loi de $N$ sous $\\mathbf{P}_\\theta$ est la loi de Poisson de paramètre $\\tilde \\lambda = \\lambda e^\\theta$. Ainsi on a la représentation \n", "\\begin{equation*}\n", " \\mathbf{P}\\biggl[ \\sum_{i=1}^N X_i > K \\biggr] \n", " = \\mathbf{E}_{\\mathbf{P}_\\theta} \\biggl[ \n", " \\mathbf{1}_{\\sum_{i=1}^N X_i > K } \\exp\\big(-\\theta N+\\psi(\\theta)\\big)\n", " \\bigg]\n", " \\quad \\text{avec $N \\sim \\mathcal{P}(\\tilde \\lambda)$ sous $\\mathbf{P}_\\theta$}.\n", "\\end{equation*}\n", "Il est d'usage pour la loi de Poisson d'écrire la variable $L_\\theta$ à partir de $\\lambda$ et $\\tilde \\lambda$ (la valeur du paramètre de la loi de Poisson sous la nouvelle probabilité) _i.e._ \n", "\\begin{equation*}\n", " L_\\theta = \\exp\\big(-\\theta N + \\lambda(e^\\theta - 1)\\big) \n", " = \\Bigl(\\frac{\\lambda}{\\tilde \\lambda}\\Big)^N \\exp\\big(\\tilde \\lambda - \\lambda\\big).\n", "\\end{equation*}" ] }, { "cell_type": "markdown", "id": "5ccb3698-b3c9-4b5d-a468-acbf3a9c41b7", "metadata": { "tags": [ "question" ] }, "source": [ "### Question: simulation sous $\\mathbf{P}_\\theta$ \n", "\n", "La loi de $N$ sous $\\mathbf{P}_\\theta$ est la loi de Poisson de paramètre $\\tilde \\lambda=\\lambda e^\\theta$ et la suite $(X_i)_{i \\ge 1}$ est indépendante de $N$ donc de $L_\\theta$ et n'est donc pas impactée par le changement de probabilité: la loi des $(X_i)_{i \\ge 1}$ est inchangée. \n", "Ecrire une fonction `simu_S_tilde` inspirée de `simu_S` qui prend un paramètre supplémentaire $\\theta$ et qui renvoie un échantillon de $\\sum_{i=1}^{N} X_i$ sous $\\mathbf{P}_\\theta$." ] }, { "cell_type": "code", "execution_count": null, "id": "4d66df95-b6d8-48a9-bbd8-a8aa955b200a", "metadata": { "tags": [ "aremplir" ] }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "e33c3558-5537-4bac-bf71-0aab45a937b5", "metadata": { "tags": [ "question" ] }, "source": [ "### Question: Monte Carlo sous $\\mathbf{P}_\\theta$\n", "\n", "Comparer pour différentes valeurs de $K$, avec $\\theta = 0.7$, l'estimateur de Monte Carlo basé sur la représentation \n", "\\begin{equation*}\n", " \\mathbf{P}\\biggl[ \\sum_{i=1}^N X_i > K \\biggr] \n", " = \\mathbf{E}\\biggl[\\mathbf{1}_{\\sum_{i=1}^{\\tilde N} X_i > K } \n", " \\bigl(\\frac{\\lambda}{\\tilde \\lambda} \\bigr)^{\\tilde N} \n", " \\exp(\\tilde \\lambda - \\lambda)\n", " \\biggr] \\quad \\text{avec $\\tilde N \\sim \\mathcal{P}(\\tilde \\lambda)$}.\n", "\\end{equation*}\n", "Pour $K = 22$ le ratio de variance est de l'ordre de 16-17.\n", "\n", "Que se passe-t-il si le paramètre $\\theta$ est mal choisi? (prendre par exemple $\\theta = 1.2$ puis $\\theta = 1.5$, et $\\theta = -0.1$...)" ] }, { "cell_type": "code", "execution_count": null, "id": "c6317f53-019c-4410-b001-9fbd797adc41", "metadata": { "tags": [ "aremplir" ] }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.10" }, "toc-autonumbering": true }, "nbformat": 4, "nbformat_minor": 5 }