Algoritmo on-line para desvio médio absoluto e grande conjunto de dados

Eu tenho um pequeno problema que está me deixando em pânico. Eu tenho que escrever o procedimento para um processo de aquisição on-line de uma série temporal multivariada. A cada intervalo de tempo (por exemplo, 1 segundo), recebo uma nova amostra, que é basicamente um vetor de ponto flutuante do tamanho N. A operação que preciso fazer é um pouco complicada:

1. Para cada nova amostra, calculo os percentuais dessa amostra (normalizando o vetor para que os elementos sejam somados a 1).

2. Calculo o vetor percentual médio da mesma maneira, mas usando os valores passados.

3. Para cada valor passado, calculo o desvio absoluto do vetor percentual relacionado a essa amostra com o vetor percentual médio global calculado na etapa 2. Dessa forma, o desvio absoluto é sempre um número entre 0 (quando o vetor é igual à média vector) e 2 (quando é totalmente diferente).

4. Usando a média dos desvios para todas as amostras anteriores, calculo o desvio absoluto médio, que é novamente um número entre 0 e 2.

5. Utilizo o desvio médio absoluto para detectar se uma nova amostra é compatível com as outras amostras (comparando seu desvio absoluto com o desvio absoluto médio de todo o conjunto calculado na etapa 4).

Como toda vez que uma nova amostra é coletada, a média global é alterada (e o desvio absoluto médio também é alterado), existe uma maneira de calcular esse valor sem verificar o conjunto de dados inteiro várias vezes? (uma vez para calcular os percentuais médios globais e uma vez para coletar os desvios absolutos). Ok, eu sei que é absolutamente fácil calcular as médias globais sem varrer todo o conjunto, já que eu só tenho que usar um vetor temporário para armazenar a soma de cada dimensão, mas e o desvio absoluto médio? Seu cálculo inclui o abs()operador, então eu preciso acessar todos os dados passados!

Obrigado pela ajuda.

Se você pode aceitar alguma imprecisão, este problema pode ser resolvido facilmente por binning contagem. Ou seja, escolha um número grande (digamos M = 1000 ) e inicialize alguns compartimentos inteiros B i , j para i = 1 … M j = 1 … N N k B i , j j ( i - 1 ) / M i /$M$$M = 1000$$B_{i,j}$$i = 1\ldots M$ e , em que é o tamanho do vetor como zero. Então, quando você observar a ésima observação de um vetor percentual, aumente se o ésimo elemento desse vetor estiver entre e , fazendo um loop sobre o N$j = 1\ldots N$$N$$k$$B_{i,j}$$j$$(i-1)/M$$i/M$$N$elementos do vetor. (Suponho que seus vetores de entrada não sejam negativos, de modo que, quando você calcula seus 'percentuais', os vetores estão no intervalo .)$[0,1]$

A qualquer momento, é possível estimar o vetor médio dos compartimentos e o desvio médio absoluto. Após observar tais vetores, o j- ésimo elemento da média é estimado por ˉ X j = $K$$j$e oj-ésimo elemento do desvio médio absoluto é estimado por

editar : este é um caso específico de uma abordagem mais geral em que você está construindo uma estimativa de densidade empírica. Isso pode ser feito com polinômios, splines, etc., mas a abordagem de binning é a mais fácil de descrever e implementar.

Eu usei a seguinte abordagem no passado para calcular o desvio de absolvição moderadamente eficiente (observe que essa é uma abordagem de programadores, não um estatístico, então, indubitavelmente, pode haver truques inteligentes como o shabbychef que podem ser mais eficientes).

AVISO: Este não é um algoritmo online. Isso requer O(n)memória. Além disso, possui o pior desempenho possível O(n), para conjuntos de dados como [1, -2, 4, -8, 16, -32, ...](ou seja, o mesmo que o recálculo completo). [1]

No entanto, como ainda funciona bem em muitos casos de uso, pode valer a pena postar aqui. Por exemplo, para calcular o desvio absoluto de 10000 números aleatórios entre -100 e 100 à medida que cada item chega, meu algoritmo leva menos de um segundo, enquanto o recálculo completo leva mais de 17 segundos (na minha máquina, variará por máquina e de acordo com os dados de entrada). No entanto, é necessário manter o vetor inteiro na memória, o que pode ser uma restrição para alguns usos. O esboço do algoritmo é o seguinte:

1. Em vez de ter um único vetor para armazenar medições passadas, use três filas de prioridade classificadas (algo como um heap mínimo / máximo). Essas três listas dividem a entrada em três: itens maiores que a média, itens menores que a média e itens iguais à média.
2. (Quase) sempre que você adiciona um item, a média muda, por isso precisamos reparticionar. O crucial é a natureza classificada das partições, o que significa que, em vez de digitalizar todos os itens da lista para reparação, precisamos apenas ler os itens que estamos movendo. Embora, no pior caso, isso ainda exija O(n)operações de movimentação, para muitos casos de uso, não é o caso.
3. Usando alguma contabilidade inteligente, podemos garantir que o desvio seja calculado corretamente o tempo todo, ao reparticionar e adicionar novos itens.

Alguns exemplos de código, em python, estão abaixo. Observe que apenas permite que itens sejam adicionados à lista, não removidos. Isso poderia ser facilmente adicionado, mas no momento em que escrevi isso não era necessário. Em vez de implementar as filas de prioridade, usei a lista classificada do excelente pacote de blist de Daniel Stutzbach , que usa internamente as árvores B + Tree .

Considere este código licenciado sob a licença MIT . Não foi significativamente otimizado ou polido, mas funcionou para mim no passado. Novas versões estarão disponíveis aqui . Deixe-me saber se você tiver alguma dúvida ou encontrar algum erro.

from blist import sortedlist
import operator

class deviance_list:
    def __init__(self):
        self.mean =  0.0
        self._old_mean = 0.0
        self._sum =  0L
        self._n =  0  #n items
        # items greater than the mean
        self._toplist =  sortedlist()
        # items less than the mean
        self._bottomlist = sortedlist(key = operator.neg)
        # Since all items in the "eq list" have the same value (self.mean) we don't need
        # to maintain an eq list, only a count
        self._eqlistlen = 0

        self._top_deviance =  0
        self._bottom_deviance =  0

    @property
    def absolute_deviance(self):
        return self._top_deviance + self._bottom_deviance

    def append(self,  n):
        # Update summary stats
        self._sum += n
        self._n +=  1
        self._old_mean =  self.mean
        self.mean =  self._sum /  float(self._n)

        # Move existing things around
        going_up = self.mean > self._old_mean
        self._rebalance(going_up)

        # Add new item to appropriate list
        if n >  self.mean:
            self._toplist.add(n)
            self._top_deviance +=  n -  self.mean
        elif n == self.mean: 
            self._eqlistlen += 1
        else:
            self._bottomlist.add(n)
            self._bottom_deviance += self.mean -  n


    def _move_eqs(self,  going_up):
        if going_up:
            self._bottomlist.update([self._old_mean] *  self._eqlistlen)
            self._bottom_deviance += (self.mean - self._old_mean) * self._eqlistlen
            self._eqlistlen = 0
        else:
            self._toplist.update([self._old_mean] *  self._eqlistlen)
            self._top_deviance += (self._old_mean - self.mean) * self._eqlistlen
            self._eqlistlen = 0


    def _rebalance(self, going_up):
        move_count,  eq_move_count = 0, 0
        if going_up:
            # increase the bottom deviance of the items already in the bottomlist
            if self.mean !=  self._old_mean:
                self._bottom_deviance += len(self._bottomlist) *  (self.mean -  self._old_mean)
                self._move_eqs(going_up)


            # transfer items from top to bottom (or eq) list, and change the deviances
            for n in iter(self._toplist):
                if n < self.mean:
                    self._top_deviance -= n -  self._old_mean
                    self._bottom_deviance += (self.mean -  n)
                    # we increment movecount and move them after the list
                    # has finished iterating so we don't modify the list during iteration
                    move_count +=  1
                elif n == self.mean:
                    self._top_deviance -= n -  self._old_mean
                    self._eqlistlen += 1
                    eq_move_count +=  1
                else:
                    break
            for _ in xrange(0,  move_count):
                self._bottomlist.add(self._toplist.pop(0))
            for _ in xrange(0,  eq_move_count):
                self._toplist.pop(0)

            # decrease the top deviance of the items remain in the toplist
            self._top_deviance -= len(self._toplist) *  (self.mean -  self._old_mean)
        else:
            if self.mean !=  self._old_mean:
                self._top_deviance += len(self._toplist) *  (self._old_mean -  self.mean)
                self._move_eqs(going_up)
            for n in iter(self._bottomlist): 
                if n > self.mean:
                    self._bottom_deviance -= self._old_mean -  n
                    self._top_deviance += n -  self.mean
                    move_count += 1
                elif n == self.mean:
                    self._bottom_deviance -= self._old_mean -  n
                    self._eqlistlen += 1
                    eq_move_count +=  1
                else:
                    break
            for _ in xrange(0,  move_count):
                    self._toplist.add(self._bottomlist.pop(0))
            for _ in xrange(0,  eq_move_count):
                self._bottomlist.pop(0)

            # decrease the bottom deviance of the items remain in the bottomlist
            self._bottom_deviance -= len(self._bottomlist) *  (self._old_mean -  self.mean)


if __name__ ==  "__main__":
    import random
    dv =  deviance_list()
    # Test against some random data,  and calculate result manually (nb. slowly) to ensure correctness
    rands = [random.randint(-100,  100) for _ in range(0,  1000)]
    ns = []
    for n in rands: 
        dv.append(n)
        ns.append(n)
        print("added:%4d,  mean:%3.2f,  oldmean:%3.2f,  mean ad:%3.2f" %
              (n, dv.mean,  dv._old_mean,  dv.absolute_deviance / dv.mean))
        assert sum(ns) == dv._sum,  "Sums not equal!"
        assert len(ns) == dv._n,  "Counts not equal!"
        m = sum(ns) / float(len(ns))
        assert m == dv.mean,  "Means not equal!"
        real_abs_dev = sum([abs(m - x) for x in ns])
        # Due to floating point imprecision, we check if the difference between the
        # two ways of calculating the asb. dev. is small rather than checking equality
        assert abs(real_abs_dev - dv.absolute_deviance) < 0.01, (
            "Absolute deviances not equal. Real:%.2f,  calc:%.2f" %  (real_abs_dev,  dv.absolute_deviance))

[1] Se os sintomas persistirem, consulte seu médico.

$X$$X$$X$$s$$s \sqrt{2 / \pi}$

MAD (x) é apenas dois cálculos medianos simultâneos, cada um dos quais pode ser disponibilizado online através do algoritmo binmediano .

Você pode encontrar o documento associado, bem como os códigos C e FORTRAN on-line aqui .

(este é apenas o uso de um truque inteligente em cima do truque inteligente de Shabbychef, para economizar memória).

Termo aditivo:

Existem vários métodos antigos de múltiplas passagens para calcular quantis. Uma abordagem popular é manter / atualizar um reservatório de tamanho determinístico de observações selecionadas aleatoriamente a partir do fluxo e calcular quantis recursivamente (veja esta revisão) neste reservatório. Essa abordagem (e relacionada) é substituída pela proposta acima.

A seguir, é apresentada uma aproximação imprecisa, embora a imprecisão dependa da distribuição dos dados de entrada. É um algoritmo online, mas apenas aproxima o desvio absoluto. Ele é baseado em um algoritmo bem conhecido para o cálculo da variação on-line, descrito por Welford na década de 1960. Seu algoritmo, traduzido para R, se parece com:

M2 <- 0
mean <- 0
n <- 0

var.online <- function(x){
    n <<- n + 1
    diff <- x - mean
    mean <<- mean + diff / n
    M2 <<- M2 + diff * (x - mean)
    variance <- M2 / (n - 1)
    return(variance)
}

Ele executa de maneira muito semelhante à função de variação interna de R:

set.seed(2099)
n.testitems <- 1000
n.tests <- 100
differences <- rep(NA, n.tests)
for (i in 1:n.tests){
        # Reset counters
        M2 <- 0
        mean <- 0
        n <- 0

        xs <- rnorm(n.testitems)
        for (j in 1:n.testitems){
                v <- var.online(xs[j])
        }

        differences[i] <- abs(v - var(xs))

}
summary(differences)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.000e+00 2.220e-16 4.996e-16 6.595e-16 9.992e-16 1.887e-15 

Modificar o algoritmo para calcular o desvio absoluto simplesmente envolve uma sqrtchamada adicional . No entanto, sqrtapresenta imprecisões que são refletidas no resultado:

absolute.deviance.online <- function(x){
    n <<- n + 1
    diff <- x - mean
    mean <<- mean + diff / n
    a.dev <<- a.dev + sqrt(diff * (x - mean))
    return(a.dev)
}

Os erros, calculados como acima, são muito maiores do que no cálculo da variação:

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.005126 0.364600 0.808000 0.958800 1.360000 3.312000 

No entanto, dependendo do seu caso de uso, essa magnitude do erro pode ser aceitável.