Simuler une loi multinomiale

Links: notebook, html ., PDF, python, slides ., GitHub

%matplotlib inline
Populating the interactive namespace from numpy and matplotlib
from jyquickhelper import add_notebook_menu
add_notebook_menu()
Plan
run previous cell, wait for 2 seconds

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.

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)
<Container object of 5 artists>
../_images/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.

import numpy.random as rnd
draw = rnd.multinomial(1000, poids)
draw / sum(draw)
array([ 0.211,  0.148,  0.15 ,  0.089,  0.402])

Pour avoir 1000 tirages plutôt que l’aggrégation des 1000 tirages :

draw = rnd.multinomial(1, poids, 1000)
draw
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.

import numpy
cum = numpy.cumsum( poids )   # voir http://docs.scipy.org/doc/numpy/reference/generated/numpy.cumsum.html
print(cum)
plt.bar( valeur, cum)
[ 0.2   0.35  0.5   0.6   1.  ]
<Container object of 5 artists>
../_images/code_multinomial_10_2.png

Cette fonction de répartition (x,F(x)) est croissante. On définit les cinq intervalles : A_i=]F(i),F(i+1)]. Pour simuler une loi multinomiale, il suffit de tirer un nombre aléatoire dans [0,1] puis de trouver l’intervalle A_i auquel il appartient. i est le résultat cherché. Le calcul de la distribution cumulée utilise une autre alternative : functools.reduce.

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]
[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.

import collections
c = collections.Counter(alea)
c
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.

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.clock()
    alea = [ simulation_multinomiale(poids_cumule) for i in range(0,10000) ]
    b = time.clock()
    print("  non trié", b-a)
    a = time.clock()
    alea = [ simulation_multinomiale(poids_cumule_trie) for i in range(0,10000) ]
    b = time.clock()
    print("  trié", b-a)
[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 K est le nombre de valeurs distinctes, les coûts fixes des deux méthodes sont les suivants :

  • non trié : O(K) (distribution cumulative)
  • trié : O(K + K\ln K) ((distribution cumulative + tri)

Qu’en est-il pour la fonction numpy.random.multinomial ?

poids_trie = list(sorted(poids, reverse=True))

for p in range(0,3):
    print("passage",p)
    a = time.clock()
    rnd.multinomial(10000, poids)
    b = time.clock()
    print("  non trié", b-a)
    a = time.clock()
    rnd.multinomial(10000, poids_trie)
    b = time.clock()
    print("  trié", b-a)
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 K est grand (K = numbre 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.

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.clock()
    alea = [ simulation_multinomiale(poids_cumule) for i in range(0,10000) ]
    b = time.clock()
    print("  simulation_multinomiale", b-a)
    a = time.clock()
    alea = [ vecteur [ random.randint(0,N) ]  for i in range(0,10000) ]
    b = time.clock()
    print("  bootstrap", b-a)
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.