<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 1A.1 - Simuler une loi multinomiale
On part d'une loi uniforme et on simule une loi multinomiale.
Step1:  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.
Step2:  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.
Step3:  Pour avoir 1000 tirages plutôt que l'aggrégation des 1000 tirages 
Step4:  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.
Step5:  Cette fonction de répartition $(x,F(x))$ est croissante. On définit les cinq intervalles 
Step6:  On vérifie que les probabilités d'apparition de chaque nombre sont celles attendues.
Step7:  Une première optimisation 
Step8:  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 
Step9:  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$ = 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.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
from jyquickhelper import add_notebook_menu
add_notebook_menu()
Explanation: 1A.1 - Simuler une loi multinomiale
On part d'une loi uniforme et on simule une loi multinomiale.
End of explanation
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)
Explanation: 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.
End of explanation
import numpy.random as rnd
draw = rnd.multinomial(1000, poids)
draw / sum(draw)
Explanation: 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.
End of explanation
draw = rnd.multinomial(1, poids, 1000)
draw
Explanation: Pour avoir 1000 tirages plutôt que l'aggrégation des 1000 tirages :
End of explanation
import numpy
cum = numpy.cumsum( poids )   # voir http://docs.scipy.org/doc/numpy/reference/generated/numpy.cumsum.html
print(cum)
plt.bar( valeur, cum)
Explanation: 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.
End of explanation
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]
Explanation: 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.
End of explanation
import collections
c = collections.Counter(alea)
c
Explanation: On vérifie que les probabilités d'apparition de chaque nombre sont celles attendues.
End of explanation
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)
Explanation: 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.
End of explanation
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)
Explanation: 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 ?
End of explanation
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)
Explanation: 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$ = 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.
End of explanation
<END_TASK>