Estimativa média robusta com eficiência de atualização O (1)

9

Estou procurando uma estimativa robusta da média que tem uma propriedade específica. Eu tenho um conjunto de elementos para os quais quero calcular esta estatística. Em seguida, adiciono novos elementos, um de cada vez, e para cada elemento adicional eu gostaria de recalcular a estatística (também conhecida como algoritmo online). Eu gostaria que esse cálculo de atualização fosse rápido, de preferência O (1), ou seja, não depende do tamanho da lista.

A média usual possui essa propriedade, que pode ser atualizada com eficiência, mas não é robusta para os valores extremos. Estimadores robustos típicos da média, como média inter-quartil e média aparada, não podem ser atualizados com eficiência (pois exigem a manutenção de uma lista classificada).

Eu gostaria de receber sugestões de estatísticas robustas que possam ser calculadas / atualizadas com eficiência.

Bit a bit
fonte
Por que não usar apenas um segmento inicial dos dados - como os 100 primeiros, os 1000 primeiros ou o que quer que seja - para erguer "cercas" para rastrear discrepâncias? Você não precisa atualizá-los novamente, portanto, não há necessidade de manter estruturas de dados adicionais.
whuber
@whuber Não posso garantir que a amostra inicial represente o restante dos dados. Por exemplo, a ordem na qual recebi os dados não é aleatória (imagine um cenário em que primeiro recebam valores mais altos e depois valores mais baixos).
Bitwise
11
Essa é uma observação crucial. Isso implica que você precisa tomar mais cuidado do que o habitual, porque inicialmente você obterá uma estimativa "robusta" dos valores extremos médios altos. Continuando a atualizar essa estimativa, você pode acabar jogando fora todos os valores mais baixos. Assim, você precisará de uma estrutura de dados na qual partes importantes de toda a distribuição de dados sejam registradas e atualizadas periodicamente. Confira nossos tópicos com as palavras-chave "online" e "quantil" para obter ideias. Dois desses promissores estão em stats.stackexchange.com/questions/3372 e stats.stackexchange.com/q/3377 .
whuber
Gostaria de oferecer uma recompensa, mas eu não tenho reputação suficiente
Jason S
11
Para continuar com a idéia no primeiro comentário do @ whuber, você pode manter um subconjunto aleatório de tamanho ou 1000, com amostra uniforme de todos os dados vistos até o momento. Este conjunto e as "cercas" associadas podem ser atualizadas em O (1). 1001000
Innuo 22/04

Respostas:

4

Esta solução implementa uma sugestão feita por @Innuo em um comentário à pergunta:

Você pode manter um subconjunto aleatório de tamanho 100 ou 1000 de amostra uniforme de todos os dados vistos até o momento. Este conjunto e as "cercas" associadas podem ser atualizadas em .O(1)

Depois que soubermos manter esse subconjunto, poderemos selecionar qualquer método que gostemos de estimar a média de uma população dessa amostra. Este é um método universal, sem nenhuma suposição, que funcionará com qualquer fluxo de entrada com uma precisão que pode ser prevista usando fórmulas de amostragem estatística padrão. (A precisão é inversamente proporcional à raiz quadrada do tamanho da amostra.)


Este algoritmo aceita como entrada um fluxo de dados t = 1 , 2 , , um tamanho de amostra m e gera um fluxo de amostras s ( t ), cada um dos quais representa a população X ( t ) = ( x ( 1 ) , x ( 2 ) , , x ( t ) ) . Especificamente, para 1 i x(t), t=1,2,,ms(t)X(t)=(x(1),x(2),,x(t)) , s ( i ) é uma amostra aleatória simples de tamanho m de X ( t ) (sem substituição).1its(i)mX(t)

Para que isso aconteça, basta que cada subconjunto de elementos de { 1 , 2 , , t } tenha chances iguais de serem os índices de x em s ( t ) . Isso implica a chance de que x ( i ) , 1 i < t , seja em s ( t ) igual a m / t, desde que t m .m{1,2,,t}xs(t)x(i), 1i<t,s(t)m/ttm

No começo, apenas coletamos o fluxo até que elementos sejam armazenados. Nesse ponto, existe apenas uma amostra possível; portanto, a condição de probabilidade é trivialmente satisfeita.m

O algoritmo assume quando . Suponha indutivamente que s (t=m+1 é uma amostra aleatória simples de X ( t ) para t > m . Defina provisoriamente s ( t + 1 ) = s ( t ) . Seja U ( t + 1 ) uma variável aleatória uniforme (independente de qualquer variável anterior usada para construir s ( t ) ). E ses(t)X(t)t>ms(t+1)=s(t)U(t+1)s(t) substitui um elemento escolhido aleatoriamente de s por x ( t + 1 ) . U(t+1)m/(t+1)sx(t+1) Esse é todo o procedimento!

Claramente tem probabilidade m / ( t + 1 ) de estar em s ( t + 1 ) . Além disso, pela hipótese de indução, x ( i ) tinha probabilidade m / t de estar em s ( t ) quando i t . Com probabilidade m / ( t + 1 ) × 1 / mx(t+1)m/(t+1)s(t+1)x(i)m/ts(t)itm/(t+1)×1/m= será removido de s ( t + 1 ) , de onde sua probabilidade de permanecer igual1/(t+1)s(t+1)

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

exatamente conforme necessário. Por indução, então, todas as probabilidades de inclusão do no s ( t ) estão corretas e é claro que não há correlação especial entre essas inclusões. Isso prova que o algoritmo está correto.x(i)s(t)

A eficiência do algoritmo é porque em cada estágio são computados no máximo dois números aleatórios e no máximo um elemento de uma matriz de m valores é substituído. O requisito de armazenamento é O ( m ) .O(1)mO(m)

A estrutura de dados para esse algoritmo consiste na amostra juntamente com o índice t da população X ( t ) que ele coleta. Inicialmente, pegamos s = X ( m ) e prosseguimos com o algoritmo para t = m + 1 , m + 2 , . Aqui está uma implementação para atualizar ( s , t ) com um valor x para produzir ( s , t +stX(t)s=X(m)t=m+1,m+2,.R(s,t)x . (O argumentodesempenha o papel de T m . O índice de t será mantida pelo chamador.)(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)
}

Para ilustrar e testar isso, usarei o estimador usual (não robusto) da média e compararei a média estimada em com a média real de X ( t ) (o conjunto cumulativo de dados visto em cada etapa) ) Eu escolhi um fluxo de entrada um tanto difícil que muda bastante suavemente, mas periodicamente sofre saltos dramáticos. O tamanho da amostra de m = 50 é bastante pequeno, permitindo ver flutuações da amostra nessas parcelas.s(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]))

Nesse ponto, onlineestá a sequência de estimativas médias produzida pela manutenção dessa amostra em execução de valores, enquanto é a sequência de estimativas médias produzida a partir de todos os dados disponíveis em cada momento. O gráfico mostra os dados (em cinza), (em preto) e duas aplicações independentes desse procedimento de amostragem (em cores). O contrato está dentro do erro de amostragem esperado:50actualactual

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

Figura


Para estimadores robustos da média, pesquise em nosso site termos e relacionados. Entre as possibilidades que vale a pena considerar estão as médias Winsorized e M-estimadores.

whuber
fonte
não está claro para mim como o limite de rejeição se parece nessa abordagem (por exemplo, o limite além do qual as observações são rejeitadas como outliers). Você pode adicioná-los à trama?
precisa saber é o seguinte
@ user603 O "limiar de rejeição" ou qualquer outro método robusto usado para estimar a média é irrelevante: escolha o método que você deseja estimar a média. (Nem todos os métodos robustos funcionam erigindo limites e rejeitando dados, BTW.) Isso seria feito no código da minha resposta, substituindo-o summarypor uma variante robusta.
whuber
Algo não está claro para mim neste exemplo. Os dados em cinza são "bons" ou "discrepantes". Se o anterior, parece que o ajuste está inclinado para baixo (deve ajustá-los melhor, pois a situação seria semelhante à tendência de queda do @ Bitwise que gostaríamos de seguir). Se os dados em cinza com valores mais altos de índice forem discrepantes, parece que o ajuste está inclinado para cima. Qual é o objetivo que você deseja ajustar aqui? O ajuste atual parece dividido entre esses dois cenários.
precisa saber é o seguinte
@ Death Como explicado no texto imediatamente anterior à figura, os dados em cinza são o fluxo de dados original. Sua média atual é a curva preta. As curvas coloridas são baseadas no algoritmo. Os desvios verticais das curvas coloridas em relação à curva preta são devidos à aleatoriedade na amostragem. A quantidade esperada de desvio em qualquer índice é proporcional ao desvio padrão dos valores de cinza anteriores a esse índice e inversamente proporcional à raiz quadrada do tamanho da amostra (tomada como 50 neste exemplo).
whuber
3

Você pode relacionar seu problema ao da tabela de controle recursivo. Esse gráfico de controle avaliará se uma nova observação está sob controle. Se for, essa observação é incluída na nova estimativa da média e variância (necessária para determinar os limites de controle).

Alguns antecedentes sobre gráficos de controle robustos, recursivos e univariados podem ser encontrados aqui . Um dos textos clássicos sobre controle de qualidade e gráficos de controle parece estar disponível online aqui .

Intuitivamente, usando a média, e uma variação σ 2 t - 1 como entradas, é possível determinar se uma nova observação no tempo t é um desvio em várias abordagens. Um seria declarar x t um valor atípico se estiver fora de um certo número de desvios padrão de μ t -μt1σt12txtμt1σt12), mas isso pode ter problemas se os dados não estiverem em conformidade com certas suposições distributivas. Se você quiser seguir esse caminho, suponha que você tenha determinado se um novo ponto não é estranho e gostaria de incluí-lo na sua estimativa média sem nenhuma taxa especial de esquecimento. Então você não pode fazer melhor do que:

μt=t1tμt1+1txt

Da mesma forma, você precisará atualizar a variação recursivamente:

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

No entanto, convém tentar algumas tabelas de controle mais convencionais. Outros gráficos de controle mais robustos à distribuição dos dados e ainda podem lidar com a não estacionariedade (como oμμσ2

Em relação a um gráfico como o EWMA, que esquece observações antigas e dá mais peso a novas, se você acha que seus dados são estacionários (o que significa que os parâmetros da distribuição geradora não mudam), não há necessidade de esquecer exponencialmente as observações mais antigas. Você pode definir o fator de esquecimento de acordo. No entanto, se você acha que não é estacionário, precisará selecionar um bom valor para o fator esquecimento (consulte novamente o livro didático para saber como fazer isso).

μ0 0σ0 02

Acho que uma abordagem nesse sentido levará à atualização mais rápida do seu problema.

Deathkill14
fonte
11
xt=porque(πt/106)+2t/106
@ Bitwise diz que a amostra inicial pode não representar dados futuros. Sem informações sobre o quão diferente será o restante dos dados, você basicamente não pode fazer nada. No entanto, se os dados iniciais tiverem informações sobre a não estacionariedade do processo (digamos, uma tendência de queda), novas observações poderão ser permitidas, levando em consideração o fato de que esperamos que sejam menores. No entanto, são necessárias algumas informações sobre a não estacionariedade. Você propõe um tipo patológico de não estacionariedade. Alguns métodos, por exemplo, o EWMA são ideais para um processo específico, mas geralmente são muito bons. Seu processo exigiria um trabalho mais personalizado.
precisa saber é o seguinte
(Detecto um matemático em você, porque é um movimento muito matemático descartar como "patológico" algo que você não pode lidar :-). Mas eu imploro para diferir com o seu prognóstico: métodos como os sugeridos por @Innuo podem realmente proteger contra essas "patologias" e tudo o mais que o mundo real possa lhe lançar, especialmente quando a randomização é incorporada à amostra.
whuber
Na verdade, eu concordo que não devemos descartar um problema que enfrentamos. Você poderia me vincular aos métodos discutidos pelo @Innuo (não consigo encontrá-los neste post - eles estavam nos links que você forneceu acima e eu senti falta deles?). Obrigado.
usar o seguinte comando
O(1 1)