Algoritmos para calcular variância - Algorithms for calculating variance


Da Wikipédia, a enciclopédia livre

Algoritmos para calcular variância desempenhar um papel importante na estatística computacional . Uma dificuldade fundamental na concepção de bons algoritmos para este problema é que as fórmulas para a variância pode envolver somas de quadrados, que pode levar à instabilidade numérica , bem como para estouro aritmético ao lidar com grandes valores.

algoritmo naïve

A fórmula para o cálculo da variância de um inteira população de tamanho N é:

Usando correcção de Bessel para calcular um imparcial estimativa da variância da população de um finito amostra de n observações, a fórmula é:

Portanto, um algoritmo ingênuo para calcular a variância estimada é dada pelo seguinte:

  • Vamos n ← 0, Sum ← 0, SUMSQ ← 0
  • Para cada dado x :
    • nn + 1
    • Soma ← soma + x
    • SUMSQ ← SUMSQ + x × x
  • Var = (SUMSQ - (Soma × Soma) / n) / (n - 1)

Este algoritmo pode ser facilmente adaptado para calcular a variância de uma população finita: basta dividir por N , em vez de n  - 1 na última linha.

Porque SUMSQ e (Soma × Soma) / n pode ser um número muito semelhantes, cancelamento pode levar à precisão do resultado a ser muito menor que a precisão inerente da aritmética de ponto flutuante usado para executar o cálculo. Assim, este algoritmo não deve ser usado na prática, e vários alternativo, numericamente estável, algoritmos têm sido propostos. Isto é particularmente ruim se o desvio padrão é pequena em relação à média. No entanto, o algoritmo pode ser melhorado mediante a adopção do método da média assumida .

Computing deslocou dados

Podemos usar uma propriedade da variância para evitar o cancelamento catastrófica nesta fórmula, ou seja, a variação é invariante com respeito a mudanças em um parâmetro de localização

com qualquer constante, o que leva a nova fórmula

quanto mais próximo é o valor médio mais preciso será o resultado será, mas apenas escolher um valor dentro das amostras variam vai garantir a estabilidade desejada. Se os valores são pequenos, então não há problemas com a soma de suas praças, pelo contrário, se eles são grandes necessariamente significa que a variação é grande também. Em qualquer caso, o segundo termo na fórmula é sempre menor do que o primeiro, por conseguinte, nenhum cancelamento pode ocorrer.

Se tomarmos apenas a primeira amostra de como o algoritmo pode ser escrito em linguagem de programação Python como

def shifted_data_variance(data):
   if len(data) < 2:
      return 0.0
   K = data[0]
   n = Ex = Ex2 = 0.0
   for x in data:
      n = n + 1
      Ex += x - K
      Ex2 += (x - K) * (x - K)
   variance = (Ex2 - (Ex * Ex)/n)/(n - 1)
   # use n instead of (n-1) if want to compute the exact variance of the given data
   # use (n-1) if data are samples of a larger population
   return variance

esta fórmula facilita assim o cálculo periódica, que pode ser expresso como

K = n = Ex = Ex2 = 0.0

def add_variable(x):
    if (n == 0):
      K = x
    n = n + 1
    Ex += x - K
    Ex2 += (x - K) * (x - K)

def remove_variable(x):
    n = n - 1
    Ex -= (x - K)
    Ex2 -= (x - K) * (x - K)

def get_meanvalue():
    return K + Ex / n

def get_variance():
    return (Ex2 - (Ex*Ex)/n) / (n-1)

algoritmo de duas passagens

Uma abordagem alternativa, usando uma fórmula diferente para a variância, primeiro calcula a média da amostra,

,

e, em seguida, calcula a soma dos quadrados das diferenças em relação à média,

,

onde s é o desvio padrão. Este é dado pela seguinte pseudocódigo:

def two_pass_variance(data):
    n = sum1 = sum2 = 0

    for x in data:
        n += 1
        sum1 += x

    mean = sum1 / n

    for x in data:
        sum2 += (x - mean)*(x - mean)

    variance = sum2 / (n - 1)
    return variance

Este algoritmo é numericamente estável, se n é pequeno. No entanto, os resultados de ambos os algoritmos simples ( "Naïve" e "Two-pass") pode depender excessivamente sobre a ordenação dos dados e pode dar maus resultados para grandes conjuntos de dados devido a erro de arredondamento repetido na acumulação da somas. Técnicas como a soma compensado pode ser usado para combater este erro para um grau.

algoritmo on-line de Welford

Muitas vezes, é útil para ser capaz de calcular a variação em uma única passagem, inspecionando cada valor apenas uma vez; por exemplo, quando os dados estão sendo coletados, sem armazenamento suficiente para manter todos os valores, ou quando os custos de acesso à memória dominar as de computação. Para um tal algoritmo on-line , uma relação de recorrência é necessária entre quantidades a partir do qual as estatísticas necessárias podem ser calculados de forma numericamente estável.

As seguintes fórmulas podem ser usadas para actualizar a média e (estimada) variação da sequência, para um elemento adicional, x n . Aqui, x n indica a média da amostra dos primeiros n amostras ( x 1 , ..., x n ), s 2 n a sua variação da amostra, e σ 2 n a sua variância da população.

Estas fórmulas sofrem de instabilidade numérica, como eles subtrair repetidamente um número pequeno a partir de um grande número de escalas que com n . A melhor quantidade para a atualização é a soma dos quadrados das diferenças em relação à média atual, , aqui denotado :

Este algoritmo foi encontrado por Welford, e tem sido exaustivamente analisados. É também comum para designar e .

Uma implementação exemplo Python para o algoritmo de Welford é dado abaixo.

# for a new value newValue, compute the new count, new mean, the new M2.
# mean accumulates the mean of the entire dataset
# M2 aggregates the squared distance from the mean
# count aggregates the number of samples seen so far
def update(existingAggregate, newValue):
    (count, mean, M2) = existingAggregate
    count += 1 
    delta = newValue - mean
    mean += delta / count
    delta2 = newValue - mean
    M2 += delta * delta2

    return (count, mean, M2)

# retrieve the mean, variance and sample variance from an aggregate
def finalize(existingAggregate):
    (count, mean, M2) = existingAggregate
    (mean, variance, sampleVariance) = (mean, M2/count, M2/(count - 1)) 
    if count < 2:
        return float('nan')
    else:
        return (mean, variance, sampleVariance)

Este algoritmo é muito menos propenso a perda de precisão devido ao cancelamento catastrófica , mas pode não ser tão eficiente por causa da operação de divisão dentro do loop. Para um algoritmo de duas passagens particularmente robusto para calcular a variância, pode-se calcular primeiro e subtrair uma estimativa da média, e, em seguida, usar esse algoritmo sobre os resíduos.

O algoritmo paralelo abaixo ilustra como mesclar vários conjuntos de estatísticas calculadas online.

algoritmo incrementais ponderada

O algoritmo pode ser alargado para segurar o peso das amostras desiguais, substituindo o contador simples n com a soma dos pesos visto até agora. Oeste (1979) sugere que este algoritmo incrementais :

def weighted_incremental_variance(dataWeightPairs):
    wSum = wSum2 = mean = S = 0

    for x, w in dataWeightPairs:  # Alternatively "for x, w in zip(data, weights):"
        wSum = wSum + w
        wSum2 = wSum2 + w*w
        meanOld = mean
        mean = meanOld + (w / wSum) * (x - meanOld)
        S = S + w * (x - meanOld) * (x - mean)

    population_variance = S / wSum
    # Bessel's correction for weighted samples
    # Frequency weights
    sample_frequency_variance = S / (wSum - 1)
    # Reliability weights
    sample_reliability_variance = S / (wSum - wSum2/wSum)

algoritmo paralelo

Chan et al. note que o algoritmo acima "On-line" é um caso especial de um algoritmo que funciona para qualquer partição da amostra em conjuntos , :

.

Isto pode ser útil quando, por exemplo, múltiplas unidades de processamento podem ser atribuídos a partes discretas da entrada.

O método de Chan para estimar a média é numericamente instável quando e ambos são grandes, porque o erro numérico na não está dimensionada para baixo da maneira que ele está no caso. Nesses casos, preferem .

def parallel_variance(avg_a, count_a, var_a, avg_b, count_b, var_b):
    delta = avg_b - avg_a
    m_a = var_a * (count_a - 1)
    m_b = var_b * (count_b - 1)
    M2 = m_a + m_b + delta ** 2 * count_a * count_b / (count_a + count_b)
    return M2 / (count_a + count_b - 1)

Exemplo

Suponha que todas as operações de ponto flutuante usar padrão IEEE 754 de dupla precisão aritmética. Considere o exemplo (4, 7, 13, 16) a partir de uma população infinita. Com base nesta amostra, a população estimada dizer é 10, ea estimativa não tendenciosa da variância da população é 30. Tanto o algoritmo ingênuo e de duas passagens algoritmo de calcular estes valores corretamente.

A seguir, considere a amostra ( 10 8  + 4 , 10 8  + 7 , 10 8  + 13 , 10 8  + 16 ), o que dá origem ao mesmo variância estimada como sendo a primeira amostra. O algoritmo de duas passagens calcula esta estimativa variância corretamente, mas o algoritmo ingênuo retorna 29,333333333333332 em vez de 30.

Embora esta perda de precisão pode ser tolerável e visto como uma pequena falha do algoritmo ingênuo, aumentando ainda mais o deslocamento faz com que o erro catastrófico. Considere a amostra ( 10 9  + 4 , 10 9  + 7 , 10 9  + 13 , 10 9  + 16 ). Mais uma vez a variância da população estimada de 30 é calculado corretamente pelo algoritmo de duas passagens, mas o algoritmo ingênuo agora calcula-o como -170,66666666666666. Este é um problema sério com o algoritmo ingênuo e é devido ao cancelamento catastrófica na subtração de dois números semelhantes na fase final do algoritmo.

estatísticas de ordem superior

Terriberry estende-se fórmulas de Chan para calcular o terceiro e quarto momentos centrais , necessários por exemplo, quando a estimativa de assimetria e curtose :

Aqui o são novamente as somas de potências das diferenças da média , dando

assimetria:
kurtosis:

Para o caso incremental (isto é, ), isto simplifica a:

Ao preservar o valor , apenas uma operação de divisão é necessária e as estatísticas de ordem superior pode assim ser calculado para pouco custo incremental.

Um exemplo do algoritmo on-line para kurtosis implementado conforme descrito é:

def online_kurtosis(data):
    n = mean = M2 = M3 = M4 = 0

    for x in data:
        n1 = n
        n = n + 1
        delta = x - mean
        delta_n = delta / n
        delta_n2 = delta_n * delta_n
        term1 = delta * delta_n * n1
        mean = mean + delta_n
        M4 = M4 + term1 * delta_n2 * (n*n - 3*n + 3) + 6 * delta_n2 * M2 - 4 * delta_n * M3
        M3 = M3 + term1 * delta_n * (n - 2) - 3 * delta_n * M2
        M2 = M2 + term1

    kurtosis = (n*M4) / (M2*M2) - 3
    return kurtosis

Pébaÿ estende-se ainda estes resultados à ordem arbitrária momentos centrais , para o incremento e os casos em pares, e subsequentemente Pébaÿ et al. para momentos ponderados e compostos. Pode-se também encontrar lá fórmulas semelhantes para covariância .

Choi e Sweetman oferecem dois métodos alternativos para calcular a assimetria e curtose, cada um dos quais pode salvar os requisitos de memória de computador substancial e tempo de CPU em determinadas aplicações. A primeira abordagem é para calcular os momentos estatísticos, separando os dados em caixas e, em seguida, computando os momentos da geometria do histograma resultante, que torna efectivamente um algoritmo de uma passagem para os momentos superiores. Uma vantagem é que os cálculos estatísticos momento pode ser realizada à precisão arbitrária tal que os cálculos podem ser ajustados para a precisão de, por exemplo, o formato de armazenamento de dados ou o hardware medição original. Um histograma relativo de uma variável aleatória pode ser construída da maneira convencional: a gama de valores possíveis é dividido em lixeiras e o número de ocorrências no interior de cada compartimento são contadas e representada de tal modo que a área de cada rectângulo é igual à porção dos valores das amostras dentro desse bin:

onde e representam a frequência e a frequência relativa no compartimento e é a área total do histograma. Após esta normalização, os momentos em bruto e momentos centrais de pode ser calculado a partir do histograma relativo:

onde o sobrescrito indica os momentos são calculados a partir do histograma. Para a largura bin constante destas duas expressões pode ser simplificada usando :

A segunda abordagem de Choi e Sweetman é uma metodologia analítica para combinar momentos estatísticos de segmentos individuais de um time-história de tal forma que os momentos totais resultantes são as do tempo histórico completo. Esta metodologia poderia ser utilizada para o cálculo paralelo de momentos estatísticos, com combinação subsequente desses momentos, ou por combinação de momentos estatísticos calculados em tempos sequenciais.

Se são conhecidos conjuntos de momentos estatísticos: para , em seguida, cada um pode ser expresso em termos do equivalente momentos matérias:

onde é geralmente considerado como sendo a duração do tempo-história, ou o número de pontos se é constante.

O benefício de expressar os momentos estatísticos em termos de é que os conjuntos podem ser combinados através da adição, e não há nenhum limite superior para o valor de .

onde o subscrito representa o tempo-história concatenadas ou combinados . Estes valores combinados da lata, em seguida, ser inversamente transformado em momentos matérias que representam o tempo-história completa concatenada

Relações conhecidas entre os momentos crus ( ) e os momentos centrais ( ) são então utilizados para calcular os momentos centrais do tempo-história concatenada. Finalmente, os momentos estatísticos da história concatenada são calculados a partir dos momentos centrais:

Covariância

Algoritmos muito semelhantes podem ser usados para calcular a covariância .

algoritmo naïve

O algoritmo ingênuo é:

Para o algoritmo acima, pode-se utilizar o seguinte código Python:

def naive_covariance(data1, data2):
    n = len(data1)
    sum12 = 0
    sum1 = sum(data1)
    sum2 = sum(data2)

    for i1, i2 in zip(data1, data2):
        sum12 += i1*i2

    covariance = (sum12 - sum1*sum2 / n) / n
    return covariance

Com estimativa da média

Quanto à variância, a covariância de duas variáveis aleatórias também é invariante ao deslocamento, por isso dado quaisquer dois valores constantes e ele pode ser escrito:

e, novamente, a escolha de um valor dentro do intervalo de valores irá estabilizar a fórmula contra o cancelamento catastrófica, bem como torná-lo mais robusto contra grandes somas. Tomando o primeiro valor de cada conjunto de dados, o algoritmo pode ser escrita como:

def shifted_data_covariance(dataX, dataY):
   n = len(dataX)
   if (n < 2):
     return 0
   kx = dataX[0]
   ky = dataY[0]
   Ex = Ey = Exy = 0
   for ix, iy in zip(dataX, dataY):
      Ex += ix - kx
      Ey += iy - ky
      Exy += (ix - kx) * (iy - ky)
   return (Exy - Ex * Ey / n) / n

Two-pass

O algoritmo em dois passos primeiro calcula os meios de amostragem, e, em seguida, a covariância:

O algoritmo de duas passagens pode ser escrita como:

def two_pass_covariance(data1, data2):
    n = len(data1)

    mean1 = sum(data1) / n
    mean2 = sum(data2) / n

    covariance = 0

    for i1, i2 in zip(data1, data2):
        a = i1 - mean1
        b = i2 - mean2
        covariance += a*b / n
    return covariance

Uma versão ligeiramente mais preciso compensado executa o algoritmo ingênuo completa sobre os resíduos. As quantias finais e deve ser zero, mas a segunda passagem compensa qualquer pequeno erro.

Conectados

Um algoritmo de uma passagem estável existe, semelhante ao algoritmo on-line para o cálculo da variância, que calcula co-momento :

A aparente assimetria na última equação é devido ao fato de que , assim que os dois termos de atualização são iguais . Mesmo uma maior precisão pode ser conseguido por em primeiro lugar computar os meios, em seguida, usando o algoritmo de uma passagem estável nos resíduos.

Assim, podemos calcular a covariância como

def online_covariance(data1, data2):
    meanx = meany = C = n = 0
    for x, y in zip(data1, data2):
        n += 1
        dx = x - meanx
        meanx += dx / n
        meany += (y - meany) / n
        C += dx * (y - meany)

    population_covar = C / n
    # Bessel's correction for sample variance
    sample_covar = C / (n - 1)

Nós também podemos fazer uma pequena modificação para calcular a covariância ponderada:

def online_weighted_covariance(data1, data2, data3):
    meanx = meany = 0
    wsum = wsum2 = 0
    C = 0
    for x, y, w in zip(data1, data2, data3):
        wsum += w
        wsum2 += w*w
        dx = x - meanx
        meanx += (w / wsum) * dx
        meany += (w / wsum) * (y - meany)
        C += w * dx * (y - meany)

    population_covar = C / wsum
    # Bessel's correction for sample variance
    # Frequency weights
    sample_frequency_covar = C / (wsum - 1)
    # Reliability weights
    sample_reliability_covar = C / (wsum - wsum2 / wsum)

Do mesmo modo, não existe uma fórmula para combinar as covariâncias de dois conjuntos que podem ser usados ​​para paralelizar o cálculo:

versão ponderada batched

Uma versão do algoritmo on-line ponderada que não agrupadas atualizado também existe: deixe denotam os pesos, podemos escrever

A covariância pode então ser calculado como

Veja também

Referências

  1. ^ Um b Bo Einarsson (1 de Agosto 2005). Precisão e confiabilidade em Computação Científica . SIAM. p. 47. ISBN  978-0-89871-584-2 . Retirado 17 de de Fevereiro de 2013 .
  2. ^ Um b T.F.Chan, GH Golub e RJ LeVeque (1983). " " Os algoritmos para calcular a variância da amostra: Análise e recomendações", o estatístico americano, 37" (PDF) : 242-247.
  3. ^ Um b Schubert, Erich; Gertz, Michael (2018/07/09). "Computação paralela numericamente estável de (co-) variância" . ACM: 10. doi : 10,1145 / 3.221.269,3223036 . ISBN  9781450365055 .
  4. ^ Higham, Nicholas (2002). Exactidão e Estabilidade de Numéricos algoritmos (2 ed) (Problema 1,10) . SIAM.
  5. ^ BP Welford (1962). "Nota sobre um método para calcular somas corrigidos de quadrados e produtos" . Technometrics 4 (3): 419-420.
  6. ^ Donald E. Knuth (1998). The Art of Computer Programming , volume 2: Seminumerical Algoritmos ., 3ª ed, p. 232. Boston: Addison-Wesley.
  7. ^ Chan, Tony F .; Golub, Gene H .; LeVeque, Randall J. (1983). Algoritmos para calcular a variância da amostra: Análise e Recomendações. O Estatístico americano 37, 242-247. https://www.jstor.org/stable/2683386
  8. ^ Ling, Robert F. (1974). Comparação de vários algoritmos para Meios Computing amostras e variações. Jornal da Associação Americana de Estatística, vol. 69, No. 348, 859-866. doi : 10,2307 / 2286154
  9. ^ http://www.johndcook.com/standard_deviation.html
  10. ^ DHD Oeste (1979). Communications of the ACM , 22, 9, 532-535: actualização das estimativas dos média e variância: um método melhorado
  11. ^ Chan, Tony F. ; Golub, Gene H. ; LeVeque, Randall J. (1979), "A actualizar Fórmulas e um algoritmo aos pares para a computação Desvios de exemplo." (PDF) , Relatório Técnico STAN-CS-79-773 , Departamento de Ciência da Computação da Universidade de Stanford.
  12. ^ Terriberry, Timothy B. (2007), computação de ordem superior Momentos online
  13. ^ Pébaÿ, Philippe (2008), "Fórmulas para Robust, uma passagem Computação Paralela da Covariâncias e ordem arbitrária-Moments estatísticos" (PDF) , Relatório Técnico SAND2008-6212 , Sandia National Laboratories
  14. ^ Pébaÿ, Philippe; Terriberry, Timothy; Kolla, Hemanth; Bennett, Janine (2016), "numericamente estável, fórmulas escaláveis para Paralela e Online Computação de alta ordem Momentos multivariada Central com pesos arbitrária" , Estatística Computacional , Springer
  15. ^ Um b Choi, Myoungkeun; Sweetman, Bert (2010), "Cálculo Eficiente dos Momentos estatísticos para monitoramento de integridade estrutural" , Journal of Monitorização Estrutural , 9 (1)

links externos