Algorithmes «en ligne» (itérateur) pour estimer statistique médiane, mode, asymétrie, kurtosis?

Existe-t-il un algorithme pour estimer la médiane, le mode, l'asymétrie et / ou le kurtosis d'un ensemble de valeurs, mais cela ne nécessite PAS de stocker toutes les valeurs en mémoire à la fois?

Je voudrais calculer les statistiques de base:

• moyenne: moyenne arithmétique
• variance: moyenne des écarts au carré de la moyenne
• écart type: racine carrée de la variance
• médiane: valeur qui sépare la plus grande moitié des nombres de la plus petite moitié
• mode: valeur la plus fréquente trouvée dans l'ensemble
• asymétrie: tl; dr
• kurtosis: tl; dr

Les formules de base pour calculer l'un de ces éléments sont l'arithmétique de l'école primaire, et je les connais. Il existe également de nombreuses bibliothèques de statistiques qui les implémentent.

Mon problème est le grand nombre (milliards) de valeurs dans les ensembles que je gère: en travaillant en Python, je ne peux pas simplement faire une liste ou un hachage avec des milliards d'éléments. Même si j'ai écrit cela en C, les tableaux de milliards d'éléments ne sont pas trop pratiques.

Les données ne sont pas triées. Il est produit au hasard, à la volée, par d'autres processus. La taille de chaque ensemble est très variable et les tailles ne seront pas connues à l'avance.

J'ai déjà compris comment gérer assez bien la moyenne et la variance, en parcourant chaque valeur de l'ensemble dans n'importe quel ordre. (En fait, dans mon cas, je les prends dans l'ordre dans lequel ils sont générés.) Voici l'algorithme que j'utilise, avec la permission de http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#On-line_algorithm :

• Initialisez trois variables: count, sum et sum_of_squares
• Pour chaque valeur:
  • Nombre d'incrément.
  • Ajoutez la valeur à additionner.
  • Ajoutez le carré de la valeur à sum_of_squares.
• Divisez la somme par le nombre, en la stockant comme moyenne de la variable.
• Divisez sum_of_squares par count, en stockant comme variable mean_of_squares.
• Moyenne carrée, stockée sous forme de square_of_mean.
• Soustrayez square_of_mean de mean_of_squares, en stockant comme variance.
• Moyenne et variance de sortie.

Cet algorithme «en ligne» a des faiblesses (par exemple, des problèmes de précision car sum_of_squares devient rapidement plus grand que la plage entière ou la précision flottante), mais il me donne essentiellement ce dont j'ai besoin, sans avoir à stocker chaque valeur dans chaque ensemble.

Mais je ne sais pas s'il existe des techniques similaires pour estimer les statistiques supplémentaires (médiane, mode, asymétrie, kurtosis). Je pourrais vivre avec un estimateur biaisé, ou même une méthode qui compromet la précision dans une certaine mesure, tant que la mémoire nécessaire pour traiter N valeurs est sensiblement inférieure à O (N).

M'indiquer une bibliothèque de statistiques existante aidera aussi, si la bibliothèque a des fonctions pour calculer une ou plusieurs de ces opérations "en ligne".

Skewness et Kurtosis

Pour les algorithmes en ligne pour Skewness et Kurtosis (dans le sens de la variance), voir dans la même page wiki ici les algorithmes parallèles pour les statistiques des moments supérieurs.

Médian

La médiane est difficile sans données triées. Si vous savez combien de points de données vous avez, en théorie, vous n'avez qu'à trier partiellement, par exemple en utilisant un algorithme de sélection . Cependant, cela n'aide pas trop avec des milliards de valeurs. Je suggérerais d'utiliser les comptages de fréquence, voir la section suivante.

Médiane et mode avec comptages de fréquence

S'il s'agit de nombres entiers, je compterais les fréquences , coupant probablement les valeurs les plus élevées et les plus basses au-delà d'une certaine valeur pour laquelle je suis sûr que ce n'est plus pertinent. Pour les flottants (ou trop d'entiers), je créerais probablement des buckets / intervalles, puis utiliserais la même approche que pour les entiers. Le mode (approximatif) et le calcul médian deviennent plus faciles, basés sur le tableau des fréquences.

Variables aléatoires normalement distribuées

S'il est distribué normalement, j'utiliserais la moyenne , la variance , l' asymétrie et l' aplatissement de l'échantillon de population comme estimateurs du maximum de vraisemblance pour un petit sous-ensemble. Les algorithmes (en ligne) pour les calculer, vous êtes déjà maintenant. Par exemple, lisez quelques centaines de milliers ou millions de points de données, jusqu'à ce que votre erreur d'estimation devienne suffisamment petite. Assurez-vous simplement de choisir au hasard dans votre ensemble (par exemple, ne pas introduire de biais en choisissant les 100'000 premières valeurs). La même approche peut également être utilisée pour estimer le mode et la médiane pour le cas normal (pour les deux, la moyenne de l'échantillon est un estimateur).

D'autres commentaires

Tous les algorithmes ci-dessus peuvent être exécutés en parallèle (y compris de nombreux algorithmes de tri et de sélection, par exemple QuickSort et QuickSelect), si cela aide.

J'ai toujours supposé (à l'exception de la section sur la distribution normale) que nous parlions de moments d'échantillonnage, de médiane et de mode, et non d'estimateurs de moments théoriques étant donné une distribution connue.

En général, l'échantillonnage des données (c'est-à-dire en regardant uniquement un sous-ensemble) devrait être assez réussi compte tenu de la quantité de données, tant que toutes les observations sont des réalisations de la même variable aléatoire (ont les mêmes distributions) et les moments, le mode et la médiane existe effectivement pour cette distribution. La dernière mise en garde n'est pas anodine. Par exemple, la moyenne (et tous les moments supérieurs) de la distribution de Cauchy n'existent pas. Dans ce cas, la moyenne de l'échantillon d'un «petit» sous-ensemble peut être massivement éloignée de la moyenne de l'échantillon de l'ensemble de l'échantillon.

J'utilise ces estimateurs de moyenne et de médiane incrémentaux / récursifs, qui utilisent tous deux le stockage constant:

mean += eta * (sample - mean)
median += eta * sgn(sample - median)

où eta est un petit paramètre de taux d'apprentissage (par exemple 0,001), et sgn () est la fonction signum qui renvoie l'un des {-1, 0, 1}. (Utilisez une constante eta si les données ne sont pas stationnaires et que vous souhaitez suivre les changements au fil du temps; sinon, pour les sources stationnaires, vous pouvez utiliser quelque chose comme eta = 1 / n pour l'estimateur moyen, où n est le nombre d'échantillons vus ainsi loin ... malheureusement, cela ne semble pas fonctionner pour l'estimateur médian.)

Ce type d'estimateur de moyenne incrémentale semble être utilisé partout, par exemple dans les règles d'apprentissage des réseaux neuronaux non supervisés, mais la version médiane semble beaucoup moins courante, malgré ses avantages (robustesse aux valeurs aberrantes). Il semble que la version médiane pourrait être utilisée pour remplacer l'estimateur moyen dans de nombreuses applications.

J'aimerais voir un estimateur de mode incrémental d'une forme similaire ...

MISE À JOUR

Je viens de modifier l'estimateur médian incrémental pour estimer des quantiles arbitraires. En général, une fonction quantile ( http://en.wikipedia.org/wiki/Quantile_function ) vous indique la valeur qui divise les données en deux fractions: p et 1-p. Ce qui suit estime cette valeur de manière incrémentielle:

quantile += eta * (sgn(sample - quantile) + 2.0 * p - 1.0)

La valeur p doit être comprise entre [0,1]. Cela déplace essentiellement la sortie symétrique de la fonction sgn () {-1,0,1} vers un côté, en partitionnant les échantillons de données en deux bacs de taille inégale (les fractions p et 1-p des données sont inférieures / supérieures à l'estimation quantile, respectivement). Notez que pour p = 0,5, cela se réduit à l'estimateur médian.

J'ai implémenté l' algorithme P-Square pour le calcul dynamique des quantiles et des histogrammes sans stocker les observations dans un module Python que j'ai écrit appelé LiveStats . Cela devrait résoudre votre problème assez efficacement. La bibliothèque prend en charge toutes les statistiques que vous mentionnez à l'exception du mode. Je n'ai pas encore trouvé de solution satisfaisante pour l'estimation de mode.

Ryan, je crains que vous ne faites pas le droit moyenne et la variance ... Cela est arrivé il y a quelques semaines ici . Et l'un des points forts de la version en ligne (qui s'appelle en fait la méthode de Welford) est le fait qu'elle soit particulièrement précise et stable, voir la discussion ici . L'un des points forts est le fait que vous n'avez pas besoin de stocker la somme totale ou la somme totale des carrés ...

Je ne peux penser à aucune approche en ligne du mode et de la médiane, qui semblent nécessiter de considérer la liste entière à la fois. Mais il se peut très bien qu'une approche similaire à celle de la variance et de la moyenne fonctionne également pour l'asymétrie et l'aplatissement ...

L'article Wikipédia cité dans la question contient les formules de calcul de l'asymétrie et du kurtosis en ligne.

Pour le mode - je crois - il n'y a aucun moyen de faire cela en ligne. Pourquoi? Supposons que toutes les valeurs de votre entrée soient différentes de la dernière qui duplique une précédente. Dans ce cas, vous devez vous souvenir de toutes les valeurs déjà vues dans l'entrée pour détecter que la dernière valeur duplique une valeur vue avant et en fait la plus fréquente.

Pour la médiane, c'est presque la même chose - jusqu'à la dernière entrée, vous ne savez pas quelle valeur deviendra la médiane si toutes les valeurs d'entrée sont différentes, car elle pourrait être avant ou après la médiane actuelle. Si vous connaissez la longueur de l'entrée, vous pouvez trouver la médiane sans stocker toutes les valeurs en mémoire, mais vous devrez quand même en stocker beaucoup (je suppose environ la moitié) car une mauvaise séquence d'entrée pourrait déplacer fortement la médiane dans le la seconde moitié peut-être une valeur de la première moitié de la médiane.

(Notez que je fais référence uniquement au calcul exact.)

Si vous avez des milliards de points de données, il est peu probable que vous ayez besoin de réponses exactes, par opposition à des réponses proches. En général, si vous avez des milliards de points de données, le processus sous-jacent qui les génère obéira probablement à une sorte de propriété statistique de stationnarité / ergodicité / mélange. Il peut également être important que vous prévoyiez que les distributions soient raisonnablement continues ou non.

Dans ces circonstances, il existe des algorithmes pour l' estimation en ligne, à faible mémoire, des quantiles (la médiane est un cas particulier de 0,5 quantile), ainsi que des modes, si vous n'avez pas besoin de réponses exactes. C'est un domaine de statistiques actif.

exemple d'estimation quantile: http://www.computer.org/portal/web/csdl/doi/10.1109/WSC.2006.323014

exemple d'estimation de mode: Bickel DR. Estimateurs robustes du mode et de l'asymétrie des données continues. Statistiques informatiques et analyse des données. 2002; 39: 153-163. doi: 10.1016 / S0167-9473 (01) 00057-3.

Ce sont des domaines actifs de la statistique informatique. Vous entrez dans des domaines où il n'y a pas de meilleur algorithme exact, mais une diversité d'entre eux (estimateurs statistiques, en vérité), qui ont des propriétés, des hypothèses et des performances différentes. Ce sont des mathématiques expérimentales. Il existe probablement des centaines à des milliers d'articles sur le sujet.

La dernière question est de savoir si vous avez vraiment besoin de l'asymétrie et de l'aplatissement par eux-mêmes, ou plus probablement d'autres paramètres qui peuvent être plus fiables pour caractériser la distribution de probabilité (en supposant que vous ayez une distribution de probabilité!). Attendez-vous un gaussien?

Avez-vous des moyens de nettoyer / prétraiter les données pour les rendre majoritairement gaussiennes? (par exemple, les montants des transactions financières sont souvent quelque peu gaussiens après avoir pris des logarithmes). Vous attendez-vous à des écarts-types finis? Vous attendez-vous à de grosses queues? Les quantités qui vous intéressent sont-elles dans les queues ou en vrac?

Tout le monde ne cesse de dire que vous ne pouvez pas utiliser le mode en ligne, mais ce n'est tout simplement pas vrai. Voici un article décrivant un algorithme pour faire exactement ce problème inventé en 1982 par Michael E. Fischer et Steven L. Salzberg de l'Université de Yale. De l'article:

L'algorithme de recherche de majorité utilise l'un de ses registres pour le stockage temporaire d'un seul élément du flux; cet élément est le candidat actuel pour l'élément majoritaire. Le deuxième registre est un compteur initialisé à 0. Pour chaque élément du flux, nous demandons à l'algorithme d'exécuter la routine suivante. Si le compteur lit 0, installez l'élément de flux actuel en tant que nouveau candidat à la majorité (en déplaçant tout autre élément qui pourrait déjà être dans le registre). Ensuite, si l'élément courant correspond au candidat majoritaire, incrémentez le compteur; sinon, décrémentez le compteur. À ce stade du cycle, si la partie du flux vue jusqu'à présent a un élément majoritaire, cet élément est dans le registre candidat et le compteur contient une valeur supérieure à 0. Et s'il n'y a pas d'élément majoritaire? Sans faire un second passage dans les données - ce qui n'est pas possible dans un environnement de flux - l'algorithme ne peut pas toujours donner une réponse sans ambiguïté dans cette circonstance. Il promet simplement d'identifier correctement l'élément majoritaire s'il y en a un.

Il peut également être étendu pour trouver le top N avec plus de mémoire, mais cela devrait le résoudre pour le mode.

Au final, si vous n'avez pas de connaissance paramétrique a priori de la distribution je pense que vous devez stocker toutes les valeurs.

Cela dit, à moins que vous ayez affaire à une sorte de situation pathologique, le remède (Rousseuw et Bassett 1990) pourrait bien être assez bon pour vos besoins.

Il s'agit très simplement de calculer la médiane des lots de médianes.

La médiane et le mode ne peuvent pas être calculés en ligne en utilisant uniquement l'espace constant disponible. Cependant, comme la médiane et le mode sont de toute façon plus "descriptifs" que "quantitatifs", vous pouvez les estimer, par exemple en échantillonnant l'ensemble de données.

Si les données sont distribuées normalement à long terme, vous pouvez simplement utiliser votre moyenne pour estimer la médiane.

Vous pouvez également estimer la médiane en utilisant la technique suivante: établir une estimation médiane M [i] pour chaque, disons, 1 000 000 entrées dans le flux de données de sorte que M [0] soit la médiane du premier million d'entrées, M [1] le médiane du deuxième million d'entrées, etc. Utilisez ensuite la médiane de M [0] ... M [k] comme estimateur médian. Cela économise bien sûr de l'espace et vous pouvez contrôler la quantité d'espace que vous souhaitez utiliser en "réglant" le paramètre 1 000 000. Cela peut également être généralisé de manière récursive.

OK mec, essayez ceci:

pour c ++:

double skew(double* v, unsigned long n){
    double sigma = pow(svar(v, n), 0.5);
    double mu = avg(v, n);

    double* t;
    t = new double[n];

    for(unsigned long i = 0; i < n; ++i){
        t[i] = pow((v[i] - mu)/sigma, 3);
    }

    double ret = avg(t, n);

    delete [] t;
    return ret;
}

double kurt(double* v, double n){
    double sigma = pow(svar(v, n), 0.5);
    double mu = avg(v, n);

    double* t;
    t = new double[n];

    for(unsigned long i = 0; i < n; ++i){
        t[i] = pow( ((v[i] - mu[i]) / sigma) , 4) - 3;
    }

    double ret = avg(t, n);

    delete [] t;
    return ret;
}

où vous dites que vous pouvez déjà calculer la variance de l'échantillon (svar) et la moyenne (moyenne), vous les dirigez vers vos fonctions pour faire cela.

Jetez également un œil à l'approximation de Pearson. sur un si grand ensemble de données, ce serait assez similaire. 3 (moyenne - médiane) / écart type vous avez la médiane comme max - min / 2

pour le mode flottant n'a pas de sens. on les collerait généralement dans des bacs d'une taille importante (comme 1/100 * (max - min)).

Ce problème a été résolu par Pebay et al:

https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf

J'aurais tendance à utiliser des seaux, qui pourraient être adaptatifs. La taille du godet doit être la précision dont vous avez besoin. Ensuite, à mesure que chaque point de données arrive, vous en ajoutez un au nombre du compartiment concerné. Celles-ci devraient vous donner des approximations simples de la médiane et de l'aplatissement, en comptant chaque compartiment comme sa valeur pondérée par son nombre.

Le seul problème pourrait être la perte de résolution en virgule flottante après des milliards d'opérations, c'est-à-dire en ajouter une ne change plus la valeur! Pour contourner ce problème, si la taille maximale du compartiment dépasse une certaine limite, vous pouvez en retirer un grand nombre.

for j in range (1,M):
    y=np.zeros(M) # build the vector y
    y[0]=y0

    #generate the white noise
    eps=npr.randn(M-1)*np.sqrt(var)

    #increment the y vector
    for k in range(1,T):
        y[k]=corr*y[k-1]+eps[k-1]

    yy[j]=y

list.append(y)