.. _reservoirsamplingrst: ===================================== Reservoir Sampling distribué - énoncé ===================================== .. only:: html **Links:** :download:`notebook `, :downloadlink:`html `, :download:`PDF `, :download:`python `, :downloadlink:`slides `, :githublink:`GitHub|_doc/notebooks/progf/reservoir_sampling.ipynb|*` `Reservoir Sampling `__ sur Map/Reduce. Cet algorithme est un algorithme d’échantillonnage en streaming. .. code:: ipython3 from jyquickhelper import add_notebook_menu add_notebook_menu() .. contents:: :local: sampling sur Hadoop / PIG ------------------------- Voir `SAMPLE `__ pour *PIG*, `takeSample `__ sur *Spark*, ou [sample](https://spark.apache.org/docs/1.4.0/api/java/org/apache/spark/sql/DataFrame.html#sample(boolean,%20double). En pratique, il difficile d’écrire un algorithme plus efficace à partir d’un langage haut niveau si ce langage haut niveau implémente lui-même la fonctionnalité. C’est une façon de reconnaître que c’est un besoin fréquent et qu’il est recommandé d’utiliser la fonction native. Reservoir Sampling ------------------ Le `reservoir sampling `__ est une façon de tirer un échantillon aléatoire de :math:`k` nombres parmi :math:`N` nombres qui ne nécessite pas de conserver ces :math:`N` nombres en mémoire. .. code:: ipython3 ensemble = [ "%d%s" % (i, chr(i%26 + 97)) for i in range(0,10000)] ensemble[:5] .. parsed-literal:: ['0a', '1b', '2c', '3d', '4e'] .. code:: ipython3 import random def reservoir_sampling(ensemble, k): N = len(ensemble) echantillon = [] for i, e in enumerate(ensemble): if len(echantillon) < k: echantillon.append(e) else: j = random.randint(0, i-1) if j < k: echantillon[j] = e return echantillon reservoir_sampling(ensemble, 10) .. parsed-literal:: ['7629l', '6747n', '1035v', '7538y', '1113v', '7565z', '2258w', '5101f', '1235n', '9965h'] Le *reservoir sampling* reproduit un tirage sans remise. L’algorithme habituel - on tire :math:`k` éléments parmi :math:`N`, la probabilité qu’un élément fasse partie de l’échantillon est de : .. math:: p = \frac{C_{N-1}^{k-1}}{C_N^k} = \frac{ \frac{N-1!}{k-1! N-k!} } { \frac{N!}{k! N-k!} } = \frac{ N-1! k! } { N! k-1!} = \frac{k}{N} Dans le cas du *reservoir sampling*, on doit s’assurer que chaque élément à la même probabilité de faire partie de l’échantillon. On procède par récurrence. Le résultat est vrai pour :math:`k=N`. On suppose que est vrai à l’ordre :math:`N-1`. La probabilité qu’un élément de l’ensemble fasse partie de l’échantillon est :math:`\frac{k}{N-1}`. A l’ordre :math:`N`, il y a une probabilité :math:`\frac{k}{N}` de remplacer un des éléments de l’échantillon. L’élément :math:`N` a donc une probabilité :math:`\frac{k}{N}` de faire partie de l’échantillon. Pour les éléments faisant déjà partie de l’échantillon, leur probabilité de rester est :math:`\frac{N-k}{N} + \frac{k}{N}\frac{k-1}{k}=\frac{N-1}{N}`. La probabilité qu’un élément présent dans l’échantillon soit dans l’échantillon à l’itération :math:`N` est de :math:`\frac{N-1}{N} \frac{k}{N-1}=\frac{k}{N}`. L’échantillon produit par le *reservoir sampling* a les mêmes propriétés qu’un échantillon produit avec un tirage sans remise. Cet algorithme est intéressent lorsqu’on veut échantillonner sur une grande base de données car il nécessite de ne garder en mémoire que l’échantillon. Reservoir Sampling Distribué ---------------------------- Distribuer l’algorithme paraît simple : une partie des données ira sur une machine qui en tirera un échantillon, l’autre machine tirera un échantillon aléatoire sur l’autre partie des données. Par extension, sur :math:`m` machines, on obtient :math:`m` échantillons de tailles :math:`k_1, ..., k_m` tirés sur des ensembles de tailles :math:`N_1, ..., N_m`. Il faut s’assurer que la version distribuée produit un échantillon avec les mêmes propriétés. exercice 1 : combinaison ~~~~~~~~~~~~~~~~~~~~~~~~ Comment recombiner ces échantillons aléatoires ? Comment choisir les :math:`k_i` sachant qu’on doit tirer un échantillon de taille :math:`k \leqslant \sum_{i=1}^m k_i` parmi :math:`N=\sum_{i=1}^m N_i` ? exercice 2 : script PIG, Spark ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Ecrire le script PIG correspondant à l’algorithme. Petit problème théorique ------------------------ Un `Générateur congruentiel linéaire `__ imite seulement le hasard, ce n’est pas vraiment le hasard. La suite dépend d’une graine ou *seed*. Si deux machines partent de la même graine, elle produiront la même suite. Qu’en est-il du hasard ? La solution de l’exercice précédent est-elle correcte d’un point de statistique ? Reservoir Sampling pondéré -------------------------- On serait tenter d’écrire… On suppose qu’on dispose un ensemble de points :math:`(X_i)` pondéré par :math:`\omega_i \in \mathbb{R}`. On veut tirer un échantillon de :math:`k` éléments. On procède de la même manière que pour le *reservoir sampling* non pondéré. A l’étape :math:`n-1`, on dispose d’un échantillon :math:`E_{n-1}=(X_{i_1}, ..., X_{i_k})`. On tire un nombre :math:`u` selon une loi uniforme :math:`[0,1]`. On ajoute l’élément :math:`X_n` si : .. math:: u \leqslant \frac{k\omega_n}{\sum_{l=1}^k w_{i_l}} On veut montrer que pour tout :math:`j \leqslant n` : .. math:: \mathbb{P}(X_j \in E_n) = \frac{k\omega_j}{\sum_{l=1}^n w_l} On procède par récurrence en supposant cela vrai à l’ordre :math:`n-1`, donc :math:`\mathbb{P}(X_j \in E_{n-1}) = \frac{k\omega_j}{\sum_{l=1}^k w_j}`. On vérifie que cela est vrai pour :math:`n = k`. Mais a priori, il n’y a aucune raison que cela soit vrai. **Il faut faire autrement.** On utilise l’algorithme `A-Res `__. Soit :math:`E_{n-1}=(X_{i_1}, ..., X_{i_k})` pondéré par :math:`(r_1, ..., r_k)` qui sont différents des :math:`(\omega_{i_1}, ..., \omega_{i_k})`. On considère l’élément :math:`X_n` pondéré par :math:`\omega_n`. On l’ajoute à :math:`E_{n-1}` si : .. math:: \min_{1 \leqslant l \leqslant k}r_l \leqslant u^{\frac{1}{\omega_n}} Si la condition est vérifié, on enlève l’observation :math:`X_{l^*}` associé au minimum :math:`r_{l^*} =\min_{1 \leqslant l \leqslant k}r_l` et on le remplace par le point :math:`X_n` pondéré par :math:`u^{\frac{1}{\omega_n}}`. Dans le cas uniforme, :math:`\omega_n=1`. Si on considère tous les nombres aléatoires :math:`u_i`, et qu’on les trie par ordre croissant :math:`(u_{\sigma(1)}, ...,u_{\sigma(n)})`. Dans ce cas, :math:`\min_{1 \leqslant l \leqslant k}r_l = u_{\sigma(n-k+1)}`. Comme les nombres :math:`u_i` sont tirés selon une loi uniforme, :math:`\mathbb{P}(\min_{1 \leqslant l \leqslant k}r_l \leqslant u) = \frac{k}{n}`. On retrouve l’algorithme précédent. Cette probabilité ne change pas lorsque :math:`\omega_n=C \; \forall n` car les nombres :math:`u_i` sont identiquement distribués. On s’intéresse maintenant au cas où les poids ne sont pas identiques. On associe à chaque élément :math:`X_i` de poids :math:`\omega_i` un nombre :math:`r_i=u^{\frac{1}{\omega_i}}` où :math:`u` un nombre aléatoire appartient tiré selon une loi uniforme. On trie ces :math:`r_i` par ordre croissant : :math:`(r_{\sigma(1)}, ...,r_{\sigma(n)})`. A l’itération :math:`n`, l’échantillon est :math:`X_{\sigma(n-k+1)}, ..., X_{\sigma(n)}`. La fonction de répartition d’une variable aléatoire est définie par :math:`F(x) = \mathbb{P}(X \leqslant x)`. On suppose que :math:`X\geqslant 0`, la loi de :math:`X^\alpha` est :math:`F(X^\alpha) = \mathbb{P}(X^\alpha \leqslant x) = \mathbb{P}(X \leqslant x^\frac{1}{\alpha}) = F(x^\frac{1}{\alpha})`.