.. _codemultinomialrst: =================================== 1A.1 - Simuler une loi multinomiale =================================== .. only:: html **Links:** :download:`notebook `, :downloadlink:`html `, :download:`python `, :downloadlink:`slides `, :githublink:`GitHub|_doc/notebooks/1a/code_multinomial.ipynb|*` On part d’une loi uniforme et on simule une loi multinomiale. .. code:: ipython3 %matplotlib inline .. parsed-literal:: Populating the interactive namespace from numpy and matplotlib .. code:: ipython3 from jyquickhelper import add_notebook_menu add_notebook_menu() .. contents:: :local: Une variable qui suit une `loi multinomiale `__ est une variable à valeurs entières qui prend ses valeurs dans un ensemble fini, et chacune de ces valeurs a une probabilité différente. .. code:: ipython3 import matplotlib.pyplot as plt poids = [ 0.2, 0.15, 0.15, 0.1, 0.4 ] valeur = [ 0,1,2,3,4 ] plt.figure(figsize=(8,4)) plt.bar(valeur,poids) .. parsed-literal:: .. image:: code_multinomial_4_1.png Lorsqu’on simule une telle loi, chaque valeur a une probabilité de sortir proportionnelle à chaque poids. La fonction `numpy.random.multinomial `__ permet de calculer cela. .. code:: ipython3 import numpy.random as rnd draw = rnd.multinomial(1000, poids) draw / sum(draw) .. parsed-literal:: array([ 0.211, 0.148, 0.15 , 0.089, 0.402]) Pour avoir 1000 tirages plutôt que l’aggrégation des 1000 tirages : .. code:: ipython3 draw = rnd.multinomial(1, poids, 1000) draw .. parsed-literal:: array([[0, 0, 0, 1, 0], [0, 0, 1, 0, 0], [0, 1, 0, 0, 0], ..., [0, 0, 0, 0, 1], [0, 0, 1, 0, 0], [0, 1, 0, 0, 0]]) Algorithme de simulation ~~~~~~~~~~~~~~~~~~~~~~~~ Tout d’abord, on calcule la distribution cumulée (ou `fonction de répartition `__). Le calcule proposé utilise la fonction `numpy.cumsum `__. .. code:: ipython3 import numpy cum = numpy.cumsum( poids ) # voir http://docs.scipy.org/doc/numpy/reference/generated/numpy.cumsum.html print(cum) plt.bar( valeur, cum) .. parsed-literal:: [ 0.2 0.35 0.5 0.6 1. ] .. parsed-literal:: .. image:: code_multinomial_10_2.png Cette fonction de répartition :math:`(x,F(x))` est croissante. On définit les cinq intervalles : :math:`A_i=]F(i),F(i+1)]`. Pour simuler une loi multinomiale, il suffit de tirer un nombre aléatoire dans :math:`[0,1]` puis de trouver l’intervalle :math:`A_i` auquel il appartient. :math:`i` est le résultat cherché. Le calcul de la distribution cumulée utilise une autre alternative : `functools.reduce `__. .. code:: ipython3 import functools, random def simulation_multinomiale(poids): # calcule la fonction de répartition # voir https://docs.python.org/3.4/library/functools.html#functools.reduce def agg(x,y): x.append(y) return x cum = functools.reduce(agg, poids, []) x = random.random() s = 0 i = 0 while s <= x and i < len(cum): s += cum[i] i += 1 return i-1 alea = [ simulation_multinomiale(poids) for i in range(0,1000) ] alea[:10] .. parsed-literal:: [0, 4, 2, 4, 2, 3, 4, 4, 4, 4] On vérifie que les probabilités d’apparition de chaque nombre sont celles attendues. .. code:: ipython3 import collections c = collections.Counter(alea) c .. parsed-literal:: Counter({4: 405, 0: 185, 2: 153, 1: 150, 3: 107}) Une première optimisation : tri des poids ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ L’algorithme fait intervenir le calcul de la distribution cumulée. Lorsqu’on tire un grand nombre de variable aléatoire, il est intéressant de ne faire ce calcul qu’une seule fois puisqu’il ne change jamais. De même, il est plus intéressant de mettre les valeurs de plus grand poids en premier. La boucle de la fonction ``simulation_multinomiale`` s’arrêtera plus tôt. .. code:: ipython3 def simulation_multinomiale(pc): x = random.random() s = 0 i = 0 while s <= x and i < len(pc): s += pc[i] i += 1 return i-1 def agg(x,y): x.append(y) return x poids_cumule = functools.reduce(agg, poids, []) poids_cumule_trie = functools.reduce(agg, sorted(poids, reverse=True), []) print(poids_cumule, poids_cumule_trie) import time for p in range(0,3): print("passage",p) a = time.perf_counter() alea = [ simulation_multinomiale(poids_cumule) for i in range(0,10000) ] b = time.perf_counter() print(" non trié", b-a) a = time.perf_counter() alea = [ simulation_multinomiale(poids_cumule_trie) for i in range(0,10000) ] b = time.perf_counter() print(" trié", b-a) .. parsed-literal:: [0.2, 0.15, 0.15, 0.1, 0.4] [0.4, 0.2, 0.15, 0.15, 0.1] passage 0 non trié 0.052738262983552886 trié 0.026375759856477998 passage 1 non trié 0.04119122404517839 trié 0.01982565262665048 passage 2 non trié 0.025675719017215215 trié 0.020590266567126037 La seconde version est plus rapide.Son intérêt dépend du nombre d’observations aléatoire à tirer. En effet, si :math:`K` est le nombre de valeurs distinctes, les coûts fixes des deux méthodes sont les suivants : - non trié : :math:`O(K)` (distribution cumulative) - trié : :math:`O(K + K\ln K)` ((distribution cumulative + tri) Qu’en est-il pour la fonction `numpy.random.multinomial `__ ? .. code:: ipython3 poids_trie = list(sorted(poids, reverse=True)) for p in range(0,3): print("passage",p) a = time.perf_counter() rnd.multinomial(10000, poids) b = time.perf_counter() print(" non trié", b-a) a = time.perf_counter() rnd.multinomial(10000, poids_trie) b = time.perf_counter() print(" trié", b-a) .. parsed-literal:: passage 0 non trié 0.00011246838175793528 trié 5.6020372539933305e-05 passage 1 non trié 4.96058262342558e-05 trié 4.704000753008586e-05 passage 2 non trié 7.483637568839185e-05 trié 4.8750553105492145e-05 C’est plus rapide aussi. On voit aussi que cette façon de faire est beaucoup plus rapide que la version implémentée en Python pur. Cela vient du faire que le module `numpy `__ est optimisé pour le calcul numérique et surtout implémenté en langages `C++ `__ et `Fortran `__. Optimisation bootstrap ~~~~~~~~~~~~~~~~~~~~~~ C’est une technique inspiré du `bootstrap `__ qui est un peu moins précise que la version précédente mais qui peut suffire dans bien des cas. Elle est intéressante si on tire un grand nomrbre d’observations aléatoire de la même loi et si :math:`K` est grand (:math:`K` = nombre de valeurs distinctes). L’idée consiste à générer un premier grand vecteur de nombres aléatoires puis à tirer des nombres aléatoire dans ce vecteur. .. code:: ipython3 K = 100 poids = [ 1/K for i in range(0,K) ] poids_cumule = functools.reduce(agg, poids, []) vecteur = [ simulation_multinomiale(poids_cumule) for i in range(0,100000) ] N = len(vecteur)-1 for p in range(0,3): print("passage",p) a = time.perf_counter() alea = [ simulation_multinomiale(poids_cumule) for i in range(0,10000) ] b = time.perf_counter() print(" simulation_multinomiale", b-a) a = time.perf_counter() alea = [ vecteur [ random.randint(0,N) ] for i in range(0,10000) ] b = time.perf_counter() print(" bootstrap", b-a) .. parsed-literal:: passage 0 simulation_multinomiale 0.20075264012029947 bootstrap 0.03310761256989281 passage 1 simulation_multinomiale 0.20289253282658137 bootstrap 0.033825186502781435 passage 2 simulation_multinomiale 0.22474218868592288 bootstrap 0.04704770498210564 Cette façon de faire implique le stockage d’un grand nombre de variables aléatoires dans ``vecteur``. Ce procédé est plus rapide car tirer un nombre aléatoire entier est plus rapide que de simuler une variable de loi multinomial.