Estimation moyenne robuste avec efficacité de mise à jour O (1)

9

Je recherche une estimation robuste de la moyenne qui a une propriété spécifique. J'ai un ensemble d'éléments pour lesquels je veux calculer cette statistique. Ensuite, j'ajoute de nouveaux éléments un par un, et pour chaque élément supplémentaire, je voudrais recalculer la statistique (également connue sous le nom d'algorithme en ligne). Je voudrais que ce calcul de mise à jour soit rapide, de préférence O (1), c'est-à-dire ne dépend pas de la taille de la liste.

La moyenne habituelle a cette propriété de pouvoir être mise à jour efficacement, mais n'est pas robuste aux valeurs aberrantes. Les estimateurs robustes typiques de la moyenne, comme la moyenne inter-quartile et la moyenne ajustée, ne peuvent pas être mis à jour efficacement (car ils nécessitent de maintenir une liste triée).

J'apprécierais toute suggestion de statistiques solides pouvant être calculées / mises à jour efficacement.

Au niveau du bit
la source
Pourquoi ne pas simplement utiliser un segment initial des données - comme les 100 premiers ou les 1000 premiers ou quoi que ce soit - pour ériger des «barrières» pour filtrer les valeurs aberrantes? Vous n'avez pas besoin de les mettre à jour à nouveau, il n'est donc pas nécessaire de maintenir des structures de données supplémentaires.
whuber
@whuber Je ne peux pas garantir que l'échantillon initial représentera le reste des données. Par exemple, l'ordre dans lequel on me donne les données n'est pas aléatoire (imaginez un scénario où l'on me donne d'abord des valeurs plus élevées puis des valeurs plus basses).
Bitwise
1
C'est une observation cruciale. Cela implique que vous devez faire plus attention que d'habitude, car vous obtiendrez initialement une estimation «robuste» des valeurs aberrantes moyennes élevées. En continuant à mettre à jour cette estimation, vous pourriez finir par rejeter toutes les valeurs inférieures. Ainsi, vous aurez besoin d'une structure de données dans laquelle des parties clés de la distribution entière des données sont enregistrées et mises à jour périodiquement. Consultez nos discussions avec des mots clés "en ligne" et "quantile" pour des idées. Deux de ces promesses sont à stats.stackexchange.com/questions/3372 et stats.stackexchange.com/q/3377 .
whuber
J'offrirais une prime mais je n'ai pas assez de réputation
Jason S
1
Pour continuer avec l'idée dans le premier commentaire de @ whuber, vous pouvez conserver un sous-ensemble aléatoire uniformément échantillonné de taille ou partir de toutes les données vues jusqu'à présent. Cet ensemble et les "clôtures" associées peuvent être mis à jour en temps O (1). 1001000
Innuo

Réponses:

4

Cette solution met en œuvre une suggestion faite par @Innuo dans un commentaire à la question:

Vous pouvez conserver un sous-ensemble aléatoire échantillonné uniformément de taille 100 ou 1000 à partir de toutes les données vues jusqu'à présent. Cet ensemble et les "clôtures" associées peuvent être mis à jour en temps .O(1)

Une fois que nous savons comment maintenir ce sous-ensemble, nous pouvons sélectionner n'importe quelle méthode que nous aimons pour estimer la moyenne d'une population à partir d'un tel échantillon. Il s'agit d'une méthode universelle, ne faisant aucune hypothèse, qui fonctionnera avec n'importe quel flux d'entrée avec une précision qui peut être prédite à l'aide de formules d'échantillonnage statistique standard. (La précision est inversement proportionnelle à la racine carrée de la taille de l'échantillon.)


Cet algorithme accepte en entrée un flux de données une taille d'échantillon , et génère un flux d'échantillons chacun représente la population . Plus précisément, pour , est un simple échantillon aléatoire de taille de (sans remplacement).x(t), t=1,2,,ms(t)X(t)=(x(1),x(2),,x(t))1itm X ( t )s(i)mX(t)

Pour cela, il suffit que chaque sous-ensemble de éléments de ait des chances égales d'être les indices de dans . Cela implique la chance que soit en est égal à condition que .{ 1 , 2 , , t } x s ( t ) x ( i ) , 1 i < t , s ( t ) m / t t mm{1,2,,t}xs(t)x(i), 1i<t,s(t)m/ttm

Au début, nous collectons simplement le flux jusqu'à ce que éléments aient été stockés. À ce stade, il n'y a qu'un seul échantillon possible, de sorte que la condition de probabilité est trivialement satisfaite.m

L'algorithme prend le relais lorsque . Supposons par induction que est un simple échantillon aléatoire de pour . Définir provisoirement . Soit une variable aléatoire uniforme (indépendante de toute variable précédente utilisée pour construire ). Si alors remplacez un élément choisi au hasard de par . C'est toute la procédure!t=m+1s(t)X(t)t>ms(t+1)=s(t)U(t+1)s(t)U(t+1)m/(t+1)sx(t+1)

Clairement, a une probabilité d'être en . De plus, par l'hypothèse d'induction, avait une probabilité d'être dans lorsque . Avec une probabilité = elle aura été supprimée de , d'où sa probabilité de rester égalex(t+1)m/(t+1)s(t+1)x(i)m/ts(t)itm/(t+1)×1/m1/(t+1)s(t+1)

mt(11t+1)=mt+1,

exactement au besoin. Par induction, donc, toutes les probabilités d'inclusion de dans sont correctes et il est clair qu'il n'y a pas de corrélation particulière entre ces inclusions. Cela prouve que l'algorithme est correct.x(i)s(t)

L'efficacité de l'algorithme est car à chaque étape, au plus deux nombres aléatoires sont calculés et au plus un élément d'un tableau de valeurs est remplacé. Le besoin de stockage est de .O(1)mO(m)

La structure de données de cet algorithme se compose de l'échantillon ainsi que de l'indice de la population qu'il échantillonne. Initialement, nous prenons et procédons à l'algorithme pour Voici une implémentation à mettre à jour avec une valeur à produire . (L'argument joue le rôle de et est . L'index sera maintenu par l'appelant.)stX(t)s=X(m)t=m+1,m+2,.R(s,t)x(s,t+1)ntsample.sizemt

update <- function(s, x, n, sample.size) {
  if (length(s) < sample.size) {
    s <- c(s, x)
  } else if (runif(1) <= sample.size / n) {
    i <- sample.int(length(s), 1)
    s[i] <- x
  }
  return (s)
}

Pour illustrer et tester cela, j'utiliserai l'estimateur habituel (non robuste) de la moyenne et comparerai la moyenne estimée de à la moyenne réelle de (l'ensemble cumulatif de données vues à chaque étape ). J'ai choisi un flux d'entrée quelque peu difficile qui change assez facilement mais subit périodiquement des sauts spectaculaires. La taille de l'échantillon de est assez petite, ce qui nous permet de voir les fluctuations d'échantillonnage dans ces parcelles.X ( t ) m = 50s(t)X(t)m=50

n <- 10^3
x <- sapply(1:(7*n), function(t) cos(pi*t/n) + 2*floor((1+t)/n))
n.sample <- 50
s <- x[1:(n.sample-1)]
online <- sapply(n.sample:length(x), function(i) {
  s <<- update(s, x[i], i, n.sample)
  summary(s)})
actual <- sapply(n.sample:length(x), function(i) summary(x[1:i]))

À ce stade, onlineest la séquence d'estimations moyennes produites en maintenant cet échantillon de valeurs, tandis que la séquence d'estimations moyennes est produite à partir de toutes les données disponibles à chaque instant. Le graphique montre les données (en gris), (en noir) et deux applications indépendantes de cette procédure d'échantillonnage (en couleurs). L'accord correspond à l'erreur d'échantillonnage attendue:50actualactual

plot(x, pch=".", col="Gray")
lines(1:dim(actual)[2], actual["Mean", ])
lines(1:dim(online)[2], online["Mean", ], col="Red")

Figure


Pour des estimateurs robustes de la moyenne, veuillez rechercher sur notre site les et les termes connexes. Parmi les possibilités qui méritent d'être étudiées figurent les moyennes winsorisées et les estimateurs M.

whuber
la source
il n'est pas clair pour moi à quoi ressemble le seuil de rejet dans cette approche (par exemple, le seuil au-delà duquel les observations sont rejetées en tant que valeurs aberrantes). Pouvez-vous les ajouter à l'intrigue?
user603
@ user603 Le "seuil de rejet", ou toute autre méthode robuste utilisée pour estimer la moyenne, n'a pas d'importance: choisissez la méthode que vous souhaitez estimer. (Toutes les méthodes robustes ne fonctionnent pas en érigeant des seuils et en rejetant des données, BTW.) Cela serait fait dans le code de ma réponse en le remplaçant summarypar une variante robuste.
whuber
Quelque chose n'est pas clair pour moi dans cet exemple. Les données grises sont-elles «bonnes» ou «aberrantes»? Si l'a priori, il semble que l'ajustement est biaisé vers le bas (il devrait mieux les adapter car la situation serait similaire à la tendance à la baisse de @ Bitwise que nous aimerions suivre). Si les données grises à des valeurs d'index plus élevées sont aberrantes, il semble que l'ajustement soit biaisé vers le haut. Quelle est la cible que vous souhaitez intégrer ici? L'ajustement actuel semble déchiré entre ces deux scénarios.
Deathkill14
@Death Comme expliqué dans le texte précédant immédiatement la figure, les données grises sont le flux de données d'origine. Sa moyenne courante est la courbe noire. Les courbes colorées sont basées sur l'algorithme. Les écarts verticaux des courbes colorées par rapport à la courbe noire sont dus au caractère aléatoire de l'échantillonnage. La quantité d'écart attendue à n'importe quel indice est proportionnelle à l'écart type des valeurs de gris précédant cet indice et inversement proportionnelle à la racine carrée de la taille de l'échantillon (prise comme 50 dans cet exemple).
whuber
3

Vous pourriez penser à relier votre problème à celui de la carte de contrôle récursive. Une telle carte de contrôle évaluera si une nouvelle observation est en contrôle. Si tel est le cas, cette observation est incluse dans la nouvelle estimation de la moyenne et de la variance (nécessaire pour déterminer les limites de contrôle).

Vous trouverez ici des informations sur les cartes de contrôle robustes, récursives et univariées . L'un des textes classiques sur le contrôle de la qualité et les cartes de contrôle semble être disponible en ligne ici .

Intuitivement, en utilisant la moyenne a, et une variance comme entrées, vous pouvez déterminer si une nouvelle observation au temps est une valeur aberrante par un certain nombre d'approches. L'une serait de déclarer une valeur aberrante si elle est en dehors d'un certain nombre d'écarts-types de (étant donné , mais cela peut poser des problèmes si les données non conforme à certaines hypothèses de répartition. Si vous voulez emprunter cette voie, supposons que vous ayez déterminé si un nouveau point n'est pas une valeur aberrante et que vous souhaitez l'inclure dans votre estimation moyenne sans taux spécial d'oubli. Alors vous ne pouvez pas faire mieux que: σ 2 t - 1 t x t μ t - 1 σ 2 t - 1 )μt1σt12txtμt1σt12)

μt=t1tμt1+1txt

De même, vous devrez mettre à jour la variance récursivement:

σt2=t1tσt12+1t1(xtμt)2

Cependant, vous voudrez peut-être essayer des cartes de contrôle plus conventionnelles. D'autres cartes de contrôle qui sont plus robustes à la distribution des données et peuvent toujours gérer la non-stationnarité (comme le de votre processus qui monte lentement) sont l'EWMA ou le CUSUM sont recommandés (voir le manuel lié ci-dessus pour plus de détails sur les graphiques et leurs limites de contrôle). Ces méthodes seront généralement moins exigeantes en termes de calcul qu'une méthode robuste, car elles ont l'avantage de simplement devoir comparer une seule nouvelle observation à des informations dérivées d'observations non aberrantes. Vous pouvez affiner vos estimations du processus à long terme et utilisé dans les calculs de limite de contrôle de ces méthodes avec les formules de mise à jour indiquées ci-dessus si vous le souhaitez.μσ 2μσ2

En ce qui concerne un graphique comme l'EWMA, qui oublie les anciennes observations et donne plus de poids aux nouvelles, si vous pensez que vos données sont stationnaires (ce qui signifie que les paramètres de la distribution génératrice ne changent pas), il n'est pas nécessaire d'oublier de manière exponentielle les observations plus anciennes. Vous pouvez définir le facteur d'oubli en conséquence. Cependant, si vous pensez qu'il s'agit d'une non-stationnarité, vous devrez sélectionner une bonne valeur pour le facteur d'oubli (voir à nouveau le manuel pour savoir comment le faire).

Je dois également mentionner qu'avant de commencer la surveillance et l'ajout de nouvelles observations en ligne, vous devrez obtenir des estimations de et (les valeurs de paramètres initiales basées sur un ensemble de données d'apprentissage) qui ne sont pas influencées par les valeurs aberrantes. Si vous pensez qu'il y a des valeurs aberrantes dans vos données d'entraînement, vous pouvez payer le coût ponctuel de l'utilisation d'une méthode robuste pour les estimer.σ 2 0μ0σ02

Je pense qu'une approche dans ce sens conduira à la mise à jour la plus rapide de votre problème.

Deathkill14
la source
1
L'utilisation de cartes de contrôle est une idée intéressante. Il semble cependant qu'il pourrait être difficile de surmonter les défis décrits dans les commentaires sur la question. Dans le cas non stationnaire, si vous "oubliez" des valeurs plus anciennes, il semble possible que les estimations soient fortement biaisées. Par exemple, comment vos suggestions fonctionneraient-elles sur un flux de données donné par ? (Cela baisse très progressivement, saute et monte très progressivement, saute à nouveau, etc.)xt=cos(πt/106)+2t/106
whuber
@Bitwise indique que l'échantillon initial peut ne pas représenter des données futures. Sans informations sur la différence du reste des données, vous ne pouvez essentiellement rien faire. Cependant, si les données initiales contiennent des informations sur la non-stationnarité du processus (par exemple une tendance à la baisse), de nouvelles observations peuvent être admises en tenant compte du fait que nous nous attendons à ce qu'elles soient plus faibles. Cependant, quelques informations sur la non-stationnarité sont nécessaires. Vous proposez un type pathologique de non-stationnarité. Certaines méthodes, par exemple l'EWMA, sont optimales pour un processus spécifique mais sont généralement assez bonnes. Votre processus nécessiterait un travail plus personnalisé.
Deathkill14
(Je détecte un mathématicien en vous, car c'est une décision très mathématique de rejeter comme quelque chose de "pathologique" que vous ne pouvez pas gérer :-). Mais je vous prie de ne pas être d'accord avec votre pronostic: des méthodes comme celles suggérées par @Innuo peuvent en effet protéger contre de telles "pathologies" et tout le reste que le monde réel pourrait vous lancer, en particulier lorsque la randomisation est incorporée dans l'échantillonnage.
whuber
En fait, je suis d'accord qu'il ne faut pas écarter un problème auquel on est confronté. Pourriez-vous s'il vous plaît me lier aux méthodes discutées @Innuo (je ne les trouve pas dans cet article - étaient-elles dans les liens que vous avez fournis ci-dessus et je les ai manquées?) Je vous remercie.
Deathkill14
@Innuo a publié un bref commentaire sur stats.stackexchange.com/questions/56494/… suggérant qu'un échantillon aléatoire uniforme de toutes les données précédemment observées pourrait être conservé en . Bien que l'on ne sache pas exactement comment cela serait fait, le coupler avec un estimateur robuste de la moyenne constituerait une solution universelle, applicable à tout flux de données que ce soit. O(1)
whuber