Estimation en ligne des quartiles sans enregistrer les observations

J'ai besoin de calculer les quartiles (Q1, médiane et Q3) en temps réel sur un grand ensemble de données sans stocker les observations. J'ai d'abord essayé l'algorithme P carré (Jain / Chlamtac) mais je n'en étais pas satisfait (un peu trop d'utilisation du processeur et pas convaincu par la précision au moins sur mon jeu de données).

J'utilise maintenant l'algorithme FAME ( Feldman / Shavitt ) pour estimer la médiane à la volée et essayer de dériver l'algorithme pour calculer également Q1 et Q3:

M = Q1 = Q3 = first data value 
step =step_Q1 = step_Q3 = a small value
for each new data :
        # update median M 
        if M > data:
            M = M - step
        elif M < data:
            M = M + step
        if abs(data-M) < step:
            step = step /2

        # estimate Q1 using M
        if data < M:
            if Q1 > data:
                Q1 = Q1 - step_Q1
            elif Q1 < data:
                Q1 = Q1 + step_Q1
            if abs(data - Q1) < step_Q1:
                step_Q1 = step_Q1/2
        # estimate Q3 using M
        elif data > M:
            if Q3 > data:
                Q3 = Q3 - step_Q3
            elif Q3 < data:
                Q3 = Q3 + step_Q3
            if abs(data-Q3) < step_Q3:
                step_Q3 = step_Q3 /2

Pour reprendre, il utilise simplement la médiane M obtenue à la volée pour diviser l'ensemble de données en deux, puis réutiliser le même algorithme pour Q1 et Q3.

Cela semble fonctionner d'une manière ou d'une autre, mais je ne suis pas en mesure de démontrer (je ne suis pas mathématicien). Est-il défectueux? J'apprécierais toute suggestion ou toute autre technique éventuelle adaptée au problème.

Merci beaucoup pour votre aide !

==== EDIT =====

Pour ceux qui sont intéressés par de telles questions, après quelques semaines, j'ai finalement fini par simplement utiliser Reservoir Sampling avec un revervoir de 100 valeurs et cela a donné des résultats très satisfaisants (pour moi).

La médiane est le point auquel la moitié des observations tombent en dessous et la moitié au-dessus. De même, le 25e centile est la médiane des données entre le min et la médiane, et le 75e centile est la médiane entre la médiane et le max, donc oui, je pense que vous êtes sur une base solide en appliquant l'algorithme médian que vous utilisez en premier l'ensemble de données pour le partitionner, puis sur les deux pièces résultantes.

Mise à jour :

Cette question sur stackoverflow conduit à cet article: Raj Jain, Imrich Chlamtac: L'algorithme P² pour le calcul dynamique des quantiiles et des histogrammes sans stocker les observations. Commun. ACM 28 (10): 1076-1085 (1985) dont le résumé indique qu'il vous intéresse probablement beaucoup:

Un algorithme heuristique est proposé pour le calcul dynamique qf de la médiane et d'autres quantiles. Les estimations sont produites dynamiquement à mesure que les observations sont générées. Les observations ne sont pas stockées; par conséquent, l'algorithme a une exigence de stockage très petite et fixe quel que soit le nombre d'observations. Cela le rend idéal pour l'implémentation dans une puce quantile qui peut être utilisée dans les contrôleurs et enregistreurs industriels. L'algorithme est en outre étendu au tracé d'histogramme. La précision de l'algorithme est analysée.

Un très léger changement dans la méthode que vous avez publiée et vous pouvez calculer n'importe quel centile arbitraire, sans avoir à calculer tous les quantiles. Voici le code Python:

class RunningPercentile:
    def __init__(self, percentile=0.5, step=0.1):
        self.step = step
        self.step_up = 1.0 - percentile
        self.step_down = percentile
        self.x = None

    def push(self, observation):
        if self.x is None:
            self.x = observation
            return

        if self.x > observation:
            self.x -= self.step * self.step_up
        elif self.x < observation:
            self.x += self.step * self.step_down
        if abs(observation - self.x) < self.step:
            self.step /= 2.0

et un exemple:

import numpy as np
import matplotlib.pyplot as plt

distribution = np.random.normal
running_percentile = RunningPercentile(0.841)
observations = []
for _ in range(1000000):
    observation = distribution()
    running_percentile.push(observation)
    observations.append(observation)

plt.figure(figsize=(10, 3))
plt.hist(observations, bins=100)
plt.axvline(running_percentile.x, c='k')
plt.show()