Algoritmos para calcular a variância - Algorithms for calculating variance

Algoritmos para calcular a variância desempenham um papel importante nas estatísticas computacionais . Uma dificuldade chave no projeto de bons algoritmos para esse problema é que as fórmulas para a variância podem envolver somas de quadrados, o que pode levar à instabilidade numérica , bem como ao estouro aritmético ao lidar com valores grandes.

Algoritmo ingênuo

Uma fórmula para calcular a variância de uma população inteira de tamanho N é:

${\ displaystyle \ sigma ^ {2} = {\ overline {(x ^ {2})}} - {\ bar {x}} ^ {2} = \ displaystyle {\ frac {\ sum _ {i = 1} ^ {N} x_ {i} ^ {2} - (\ sum _ {i = 1} ^ {N} x_ {i}) ^ {2} / N} {N}}. \!}$

Usando a correção de Bessel para calcular uma estimativa imparcial da variância da população a partir de uma amostra finita de n observações, a fórmula é:

${\ displaystyle s ^ {2} = \ left ({\ frac {\ sum _ {i = 1} ^ {n} x_ {i} ^ {2}} {n}} - \ left ({\ frac {\ soma _ {i = 1} ^ {n} x_ {i}} {n}} \ direita) ^ {2} \ direita) \ cdot {\ frac {n} {n-1}}. \!}$

Portanto, um algoritmo ingênuo para calcular a variação estimada é dado pelo seguinte:

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

Como SumSq e (Sum × Sum) / n podem ser números muito semelhantes, o cancelamento pode fazer com que a precisão do resultado seja muito menor do que a precisão inerente da aritmética de ponto flutuante usada para realizar o cálculo. Portanto, este algoritmo não deve ser usado na prática, e vários algoritmos alternativos, numericamente estáveis, foram propostos. Isso é particularmente ruim se o desvio padrão for pequeno em relação à média.

Computando dados deslocados

A variância é invariável em relação às mudanças em um parâmetro de localização , uma propriedade que pode ser usada para evitar o cancelamento catastrófico nesta fórmula.

${\ displaystyle \ operatorname {Var} (XK) = \ operatorname {Var} (X).}$

com qualquer constante, o que leva à nova fórmula ${\ displaystyle K}$

${\ displaystyle \ sigma ^ {2} = \ displaystyle {\ frac {\ sum _ {i = 1} ^ {n} (x_ {i} -K) ^ {2} - (\ sum _ {i = 1} ^ {n} (x_ {i} -K)) ^ {2} / n} {n-1}}. \!}$

quanto mais próximo do valor médio, mais preciso será o resultado, mas apenas a escolha de um valor dentro da faixa de amostras garantirá a estabilidade desejada. Se os valores forem pequenos, então não há problemas com a soma dos seus quadrados, pelo contrário, se forem grandes significa necessariamente que a variância também é grande. Em qualquer caso, o segundo termo na fórmula é sempre menor do que o primeiro, portanto, nenhum cancelamento pode ocorrer. ${\ displaystyle K}$${\ displaystyle (x_ {i} -K)}$

Se apenas a primeira amostra for obtida, o algoritmo pode ser escrito na linguagem de programação Python como ${\ displaystyle K}$

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 também facilita o cálculo incremental que pode ser expresso como

K = n = Ex = Ex2 = 0.0

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

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

def get_mean():
    global K, n, Ex
    return K + Ex / n

def get_variance():
    global n, Ex, Ex2
    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,

${\ displaystyle {\ bar {x}} = {\ frac {\ sum _ {j = 1} ^ {n} x_ {j}} {n}},}$

e então calcula a soma dos quadrados das diferenças da média,

${\ displaystyle {\ text {sample variance}} = s ^ {2} = \ displaystyle {\ frac {\ sum _ {i = 1} ^ {n} (x_ {i} - {\ bar {x}}) ^ {2}} {n-1}}, \!}$

onde s é o desvio padrão. Isso é dado pelo seguinte có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 for pequeno. No entanto, os resultados de ambos os algoritmos simples ("ingênuo" e "duas passagens") podem depender excessivamente da ordem dos dados e podem dar resultados ruins para conjuntos de dados muito grandes devido a repetidos erros de arredondamento no acúmulo de somas. Técnicas como soma compensada podem ser usadas para combater esse erro até certo ponto.

Algoritmo online de Welford

Muitas vezes é útil 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 dominam os de computação. Para tal algoritmo online , uma relação de recorrência é necessária entre as quantidades a partir das quais as estatísticas necessárias podem ser calculadas de forma numericamente estável. ${\ displaystyle x_ {i}}$

As fórmulas a seguir podem ser usadas para atualizar a média e a variância (estimada) da sequência, para um elemento adicional x _n . Aqui, denota a média da amostra das primeiras n amostras , sua variância de amostra enviesada e sua variância de amostra enviesada . ${\ displaystyle {\ overline {x}} _ {n} = {\ frac {1} {n}} \ sum _ {i = 1} ^ {n} x_ {i}}$${\ displaystyle (x_ {1}, \ dots, x_ {n})}$${\ displaystyle \ sigma _ {n} ^ {2} = {\ frac {1} {n}} \ sum _ {i = 1} ^ {n} \ left (x_ {i} - {\ overline {x} } _ {n} \ right) ^ {2}}$${\ displaystyle s_ {n} ^ {2} = {\ frac {1} {n-1}} \ sum _ {i = 1} ^ {n} \ left (x_ {i} - {\ overline {x} } _ {n} \ right) ^ {2}}$

${\ displaystyle {\ bar {x}} _ {n} = {\ frac {(n-1) \, {\ bar {x}} _ {n-1} + x_ {n}} {n}} = {\ bar {x}} _ {n-1} + {\ frac {x_ {n} - {\ bar {x}} _ {n-1}} {n}} \!}$

${\ displaystyle \ sigma _ {n} ^ {2} = {\ frac {(n-1) \, \ sigma _ {n-1} ^ {2} + (x_ {n} - {\ bar {x} } _ {n-1}) (x_ {n} - {\ bar {x}} _ {n})} {n}} = \ sigma _ {n-1} ^ {2} + {\ frac {( x_ {n} - {\ bar {x}} _ {n-1}) (x_ {n} - {\ bar {x}} _ {n}) - \ sigma _ {n-1} ^ {2} } {n}}.}$

${\ displaystyle s_ {n} ^ {2} = {\ frac {n-2} {n-1}} \, s_ {n-1} ^ {2} + {\ frac {(x_ {n} - { \ bar {x}} _ {n-1}) ^ {2}} {n}} = s_ {n-1} ^ {2} + {\ frac {(x_ {n} - {\ bar {x} } _ {n-1}) ^ {2}} {n}} - {\ frac {s_ {n-1} ^ {2}} {n-1}}, \ quad n> 1}$

Essas fórmulas sofrem de instabilidade numérica, pois repetidamente subtraem um pequeno número de um grande número que escala com n . Uma quantidade melhor para atualização é a soma dos quadrados das diferenças da média atual , aqui denotada : ${\ displaystyle \ textstyle \ sum _ {i = 1} ^ {n} (x_ {i} - {\ bar {x}} _ {n}) ^ {2}}$${\ displaystyle M_ {2, n}}$

${\ displaystyle {\ begin {alinhados} M_ {2, n} & = M_ {2, n-1} + (x_ {n} - {\ bar {x}} _ {n-1}) (x_ {n } - {\ bar {x}} _ {n}) \\ [4pt] \ sigma _ {n} ^ {2} & = {\ frac {M_ {2, n}} {n}} \\ [4pt ] s_ {n} ^ {2} & = {\ frac {M_ {2, n}} {n-1}} \ end {alinhado}}}$

Este algoritmo foi encontrado por Welford e foi completamente analisado. Também é comum denotar e . ${\ displaystyle M_ {k} = {\ bar {x}} _ {k}}$${\ displaystyle S_ {k} = M_ {2, k}}$

Um exemplo de implementação Python para o algoritmo de Welford é fornecido 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
    if count < 2:
        return float("nan")
    else:
        (mean, variance, sampleVariance) = (mean, M2 / count, M2 / (count - 1))
        return (mean, variance, sampleVariance)

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

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

Algoritmo incremental ponderado

O algoritmo pode ser estendido para lidar com pesos amostrais desiguais, substituindo o contador simples n pela soma dos pesos vistos até agora. West (1979) sugere este algoritmo incremental :

def weighted_incremental_variance(data_weight_pairs):
    w_sum = w_sum2 = mean = S = 0

    for x, w in data_weight_pairs:
        w_sum = w_sum + w
        w_sum2 = w_sum2 + w * w
        mean_old = mean
        mean = mean_old + (w / w_sum) * (x - mean_old)
        S = S + w * (x - mean_old) * (x - mean)

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

Algoritmo paralelo

Chan et al. observe que o algoritmo online de Welford detalhado acima é um caso especial de um algoritmo que funciona para combinar conjuntos arbitrários e : ${\ displaystyle A}$${\ displaystyle B}$

${\ displaystyle {\ begin {alinhado} n_ {AB} & = n_ {A} + n_ {B} \\\ delta & = {\ bar {x}} _ {B} - {\ bar {x}} _ {A} \\ {\ bar {x}} _ {AB} & = {\ bar {x}} _ {A} + \ delta \ cdot {\ frac {n_ {B}} {n_ {AB}}} \\ M_ {2, AB} & = M_ {2, A} + M_ {2, B} + \ delta ^ {2} \ cdot {\ frac {n_ {A} n_ {B}} {n_ {AB} }} \\\ end {alinhado}}}$.

Isso pode ser útil quando, por exemplo, várias unidades de processamento podem ser atribuídas 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 em não é reduzido da maneira que é no caso. Nesses casos, prefira . ${\ displaystyle n_ {A} \ aprox n_ {B}}$${\ displaystyle \ delta = {\ bar {x}} _ {B} - {\ bar {x}} _ {A}}$${\ displaystyle n_ {B} = 1}$${\ displaystyle {\ bar {x}} _ {AB} = {\ frac {n_ {A} {\ bar {x}} _ {A} + n_ {B} {\ bar {x}} _ {B} } {n_ {AB}}}}$

def parallel_variance(n_a, avg_a, M2_a, n_b, avg_b, M2_b):
    n = n_a + n_b
    delta = avg_b - avg_a
    M2 = M2_a + M2_b + delta ** 2 * n_a * n_b / n
    var_ab = M2 / (n - 1)
    return var_ab

Isso pode ser generalizado para permitir a paralelização com AVX , com GPUs e clusters de computador e para covariância.

Exemplo

Suponha que todas as operações de ponto flutuante usem aritmética de precisão dupla padrão IEEE 754 . Considere a amostra (4, 7, 13, 16) de uma população infinita. Com base nessa amostra, a média da população estimada é 10 e a estimativa não enviesada da variância da população é 30. Tanto o algoritmo ingênuo quanto o algoritmo de duas passagens calculam esses valores corretamente.

Em seguida, considere a amostra ( 10 ⁸  + 4 , 10 ⁸  + 7 , 10 ⁸  + 13 , 10 ⁸  + 16 ), que dá origem à mesma variância estimada da primeira amostra. O algoritmo de duas passagens calcula essa estimativa de variância corretamente, mas o algoritmo ingênuo retorna 29.333333333333332 em vez de 30.

Embora essa perda de precisão possa ser tolerável e vista como uma falha secundária do algoritmo ingênuo, aumentar ainda mais o deslocamento torna o erro catastrófico. Considere a amostra ( 10 ⁹  + 4 , 10 ⁹  + 7 , 10 ⁹  + 13 , 10 ⁹  + 16 ). Novamente, a variação da população estimada de 30 é calculada corretamente pelo algoritmo de duas passagens, mas o algoritmo ingênuo agora a calcula como -170.66666666666666. Este é um problema sério com o algoritmo ingênuo e é devido ao cancelamento catastrófico na subtração de dois números semelhantes no estágio final do algoritmo.

Estatísticas de ordem superior

Terriberry estende as fórmulas de Chan para calcular o terceiro e o quarto momentos centrais , necessários, por exemplo, ao estimar assimetria e curtose :

${\ displaystyle {\ begin {alinhados} M_ {3, X} = M_ {3, A} + M_ {3, B} & {} + \ delta ^ {3} {\ frac {n_ {A} n_ {B } (n_ {A} -n_ {B})} {n_ {X} ^ {2}}} + 3 \ delta {\ frac {n_ {A} M_ {2, B} -n_ {B} M_ {2 , A}} {n_ {X}}} \\ [6pt] M_ {4, X} = M_ {4, A} + M_ {4, B} & {} + \ delta ^ {4} {\ frac { n_ {A} n_ {B} \ esquerda (n_ {A} ^ {2} -n_ {A} n_ {B} + n_ {B} ^ {2} \ direita)} {n_ {X} ^ {3} }} \\ [6pt] & {} + 6 \ delta ^ {2} {\ frac {n_ {A} ^ {2} M_ {2, B} + n_ {B} ^ {2} M_ {2, A }} {n_ {X} ^ {2}}} + 4 \ delta {\ frac {n_ {A} M_ {3, B} -n_ {B} M_ {3, A}} {n_ {X}}} \ end {alinhado}}}$

Aqui estão novamente as somas das potências das diferenças da média , dando ${\ displaystyle M_ {k}}$${\ displaystyle \ sum (x - {\ overline {x}}) ^ {k}}$

${\ displaystyle {\ begin {align} & {\ text {skewness}} = g_ {1} = {\ frac {{\ sqrt {n}} M_ {3}} {M_ {2} ^ {3/2} }}, \\ [4pt] & {\ text {curtose}} = g_ {2} = {\ frac {nM_ {4}} {M_ {2} ^ {2}}} - 3. \ end {alinhado} }}$

Para o caso incremental (ou seja, ), isso simplifica para: ${\ displaystyle B = \ {x \}}$

${\ displaystyle {\ begin {alinhado} \ delta & = xm \\ [5pt] m '& = m + {\ frac {\ delta} {n}} \\ [5pt] M_ {2}' & = M_ {2 } + \ delta ^ {2} {\ frac {n-1} {n}} \\ [5pt] M_ {3} '& = M_ {3} + \ delta ^ {3} {\ frac {(n- 1) (n-2)} {n ^ {2}}} - {\ frac {3 \ delta M_ {2}} {n}} \\ [5pt] M_ {4} '& = M_ {4} + {\ frac {\ delta ^ {4} (n-1) (n ^ {2} -3n + 3)} {n ^ {3}}} + {\ frac {6 \ delta ^ {2} M_ {2 }} {n ^ {2}}} - {\ frac {4 \ delta M_ {3}} {n}} \ end {alinhado}}}$

Ao preservar o valor , apenas uma operação de divisão é necessária e as estatísticas de ordem superior podem, portanto, ser calculadas com pouco custo incremental. ${\ displaystyle \ delta / n}$

Um exemplo do algoritmo online para curtose 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

    # Note, you may also calculate variance using M2, and skewness using M3
    # Caution: If all the inputs are the same, M2 will be 0, resulting in a division by 0.
    kurtosis = (n * M4) / (M2 * M2) - 3
    return kurtosis

Pébaÿ estende ainda mais esses resultados para momentos centrais de ordem arbitrária , para os casos incrementais e aos pares, e subsequentemente Pébaÿ et al. para momentos ponderados e compostos. Também podem ser encontradas 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 economizar requisitos substanciais de memória do computador e tempo de CPU em certos aplicativos. A primeira abordagem é computar os momentos estatísticos, separando os dados em caixas e, em seguida, computando os momentos da geometria do histograma resultante, que se torna efetivamente um algoritmo de uma passagem para momentos mais elevados. Um benefício é que os cálculos estatísticos do momento podem ser realizados com precisão arbitrária, de modo que os cálculos podem ser ajustados para a precisão, por exemplo, do formato de armazenamento de dados ou do hardware de medição original. Um histograma relativo de uma variável aleatória pode ser construído da maneira convencional: o intervalo de valores potenciais é dividido em caixas e o número de ocorrências dentro de cada caixa é contado e plotado de modo que a área de cada retângulo seja igual à parte dos valores de amostra dentro dessa caixa:

${\ displaystyle H (x_ {k}) = {\ frac {h (x_ {k})} {A}}}$

onde e representam a frequência e a frequência relativa no bin e é a área total do histograma. Após esta normalização, os momentos brutos e os momentos centrais de podem ser calculados a partir do histograma relativo: ${\ displaystyle h (x_ {k})}$${\ displaystyle H (x_ {k})}$${\ displaystyle x_ {k}}$${\ displaystyle A = \ sum _ {k = 1} ^ {K} h (x_ {k}) \, \ Delta x_ {k}}$${\ displaystyle n}$${\ displaystyle x (t)}$

${\ displaystyle m_ {n} ^ {(h)} = \ sum _ {k = 1} ^ {K} x_ {k} ^ {n} H (x_ {k}) \, \ Delta x_ {k} = {\ frac {1} {A}} \ sum _ {k = 1} ^ {K} x_ {k} ^ {n} h (x_ {k}) \, \ Delta x_ {k}}$

${\ displaystyle \ theta _ {n} ^ {(h)} = \ sum _ {k = 1} ^ {K} {\ Big (} x_ {k} -m_ {1} ^ {(h)} {\ Grande)} ^ {n} \, H (x_ {k}) \, \ Delta x_ {k} = {\ frac {1} {A}} \ sum _ {k = 1} ^ {K} {\ Grande (} x_ {k} -m_ {1} ^ {(h)} {\ Grande)} ^ {n} h (x_ {k}) \, \ Delta x_ {k}}$

onde o sobrescrito indica que os momentos são calculados a partir do histograma. Para largura de compartimento constante, essas duas expressões podem ser simplificadas usando : ${\ displaystyle ^ {(h)}}$${\ displaystyle \ Delta x_ {k} = \ Delta x}$${\ displaystyle I = A / \ Delta x}$

${\ displaystyle m_ {n} ^ {(h)} = {\ frac {1} {I}} \ sum _ {k = 1} ^ {K} x_ {k} ^ {n} \, h (x_ { k})}$

${\ displaystyle \ theta _ {n} ^ {(h)} = {\ frac {1} {I}} \ sum _ {k = 1} ^ {K} {\ Big (} x_ {k} -m_ { 1} ^ {(h)} {\ Grande)} ^ {n} h (x_ {k})}$

A segunda abordagem de Choi e Sweetman é uma metodologia analítica para combinar momentos estatísticos de segmentos individuais de uma história de tempo de modo que os momentos gerais resultantes sejam aqueles da história de tempo completa. Esta metodologia pode ser usada para cálculo paralelo de momentos estatísticos com subsequente combinação desses momentos, ou para combinação de momentos estatísticos computados em tempos sequenciais.

Se conjuntos de momentos estatísticos são conhecidos: para , então cada um pode ser expresso em termos de momentos brutos equivalentes : ${\ displaystyle Q}$${\ displaystyle (\ gamma _ {0, q}, \ mu _ {q}, \ sigma _ {q} ^ {2}, \ alpha _ {3, q}, \ alpha _ {4, q}) \ quad}$${\ displaystyle q = 1,2, \ ldots, Q}$${\ displaystyle \ gamma _ {n}}$${\ displaystyle n}$

${\ displaystyle \ gamma _ {n, q} = m_ {n, q} \ gamma _ {0, q} \ qquad \ quad {\ textrm {para}} \ quad n = 1,2,3,4 \ quad {\ text {and}} \ quad q = 1,2, \ dots, Q}$

onde geralmente é considerado a duração do histórico de tempo, ou o número de pontos se for constante. ${\ displaystyle \ gamma _ {0, q}}$${\ displaystyle q ^ {th}}$${\ displaystyle \ Delta t}$

A vantagem de expressar os momentos estatísticos em termos de é que os conjuntos podem ser combinados por adição e não há limite superior para o valor de . ${\ displaystyle \ gamma}$${\ displaystyle Q}$${\ displaystyle Q}$

${\ displaystyle \ gamma _ {n, c} = \ sum _ {q = 1} ^ {Q} \ gamma _ {n, q} \ quad \ quad {\ text {for}} n = 0,1,2 , 3,4}$

onde o subscrito representa o histórico de tempo concatenado ou combinado . Esses valores combinados de podem então ser inversamente transformados em momentos brutos que representam a história do tempo concatenada completa ${\ displaystyle _ {c}}$${\ displaystyle \ gamma}$${\ displaystyle \ gamma}$

${\ displaystyle m_ {n, c} = {\ frac {\ gamma _ {n, c}} {\ gamma _ {0, c}}} \ quad {\ text {for}} n = 1,2,3 , 4}$

Relacionamentos conhecidos entre os momentos brutos ( ) e os momentos centrais ( ) são então usados ​​para computar os momentos centrais da história do tempo concatenada. Finalmente, os momentos estatísticos da história concatenada são calculados a partir dos momentos centrais: ${\ displaystyle m_ {n}}$${\ displaystyle \ theta _ {n} = \ operatorname {E} [(x- \ mu) ^ {n}])}$

${\ displaystyle \ mu _ {c} = m_ {1, c} \ qquad \ sigma _ {c} ^ {2} = \ theta _ {2, c} \ qquad \ alpha _ {3, c} = {\ frac {\ theta _ {3, c}} {\ sigma _ {c} ^ {3}}} \ qquad \ alpha _ {4, c} = {\ frac {\ theta _ {4, c}} {\ sigma _ {c} ^ {4}}} - 3}$

Covariância

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

Algoritmo ingênuo

O algoritmo ingênuo é:

${\ displaystyle \ operatorname {Cov} (X, Y) = \ displaystyle {\ frac {\ sum _ {i = 1} ^ {n} x_ {i} y_ {i} - (\ sum _ {i = 1} ^ {n} x_ {i}) (\ sum _ {i = 1} ^ {n} y_ {i}) / n} {n}}. \!}$

Para o algoritmo acima, pode-se usar 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, portanto, dados quaisquer dois valores constantes e podem ser escritos: ${\ displaystyle k_ {x}}$${\ displaystyle k_ {y},}$

${\ displaystyle \ operatorname {Cov} (X, Y) = \ operatorname {Cov} (X-k_ {x}, Y-k_ {y}) = \ displaystyle {\ frac {\ sum _ {i = 1} ^ {n} (x_ {i} -k_ {x}) (y_ {i} -k_ {y}) - (\ sum _ {i = 1} ^ {n} (x_ {i} -k_ {x}) ) (\ sum _ {i = 1} ^ {n} (y_ {i} -k_ {y})) / n} {n}}. \!}$

e, novamente, a escolha de um valor dentro da faixa de valores estabilizará a fórmula contra o cancelamento catastrófico, bem como a tornará mais robusta contra grandes somas. Tomando o primeiro valor de cada conjunto de dados, o algoritmo pode ser escrito como:

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

Duas passagens

O algoritmo de duas passagens primeiro calcula as médias da amostra e, em seguida, a covariância:

${\ displaystyle {\ bar {x}} = \ displaystyle \ sum _ {i = 1} ^ {n} x_ {i} / n}$

${\ displaystyle {\ bar {y}} = \ displaystyle \ sum _ {i = 1} ^ {n} y_ {i} / n}$

${\ displaystyle \ operatorname {Cov} (X, Y) = \ displaystyle {\ frac {\ sum _ {i = 1} ^ {n} (x_ {i} - {\ bar {x}}) (y_ {i } - {\ bar {y}})} {n}}. \!}$

O algoritmo de duas passagens pode ser escrito 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 compensada um pouco mais precisa executa o algoritmo ingênuo completo nos resíduos. As somas finais e devem ser zero, mas a segunda passagem compensa qualquer pequeno erro. ${\ displaystyle \ textstyle \ sum x_ {i}}$${\ displaystyle \ textstyle \ sum y_ {i}}$

Conectados

Existe um algoritmo de passagem única estável, semelhante ao algoritmo online para calcular a variância, que calcula o co-momento : ${\ displaystyle \ textstyle C_ {n} = \ sum _ {i = 1} ^ {n} (x_ {i} - {\ bar {x}} _ {n}) (y_ {i} - {\ bar { y}} _ {n})}$

${\ displaystyle {\ begin {alignat} {2} {\ bar {x}} _ {n} & = {\ bar {x}} _ {n-1} & \, + \, & {\ frac {x_ {n} - {\ bar {x}} _ {n-1}} {n}} \\ [5pt] {\ bar {y}} _ {n} & = {\ bar {y}} _ {n -1} & \, + \, & {\ frac {y_ {n} - {\ bar {y}} _ {n-1}} {n}} \\ [5pt] C_ {n} & = C_ { n-1} & \, + \, & (x_ {n} - {\ bar {x}} _ {n}) (y_ {n} - {\ bar {y}} _ {n-1}) \ \ [5pt] & = C_ {n-1} & \, + \, & (x_ {n} - {\ bar {x}} _ {n-1}) (y_ {n} - {\ bar {y }} _ {n}) \ end {alignat}}}$

A aparente assimetria nessa última equação se deve ao fato de que , ambos os termos de atualização são iguais a . Uma precisão ainda maior pode ser alcançada calculando primeiro os meios e, em seguida, usando o algoritmo de uma passagem estável nos resíduos. ${\ displaystyle \ textstyle (x_ {n} - {\ bar {x}} _ {n}) = {\ frac {n-1} {n}} (x_ {n} - {\ bar {x}} _ {n-1})}$${\ displaystyle \ textstyle {\ frac {n-1} {n}} (x_ {n} - {\ bar {x}} _ {n-1}) (y_ {n} - {\ bar {y}} _ {n-1})}$

Assim, a covariância pode ser calculada como

${\ displaystyle {\ begin {alinhados} \ operatorname {Cov} _ {N} (X, Y) = {\ frac {C_ {N}} {N}} & = {\ frac {\ operatorname {Cov} _ { N-1} (X, Y) \ cdot (N-1) + (x_ {n} - {\ bar {x}} _ {n}) (y_ {n} - {\ bar {y}} _ { n-1})} {N}} \\ & = {\ frac {\ operatorname {Cov} _ {N-1} (X, Y) \ cdot (N-1) + (x_ {n} - {\ bar {x}} _ {n-1}) (y_ {n} - {\ bar {y}} _ {n})} {N}} \\ & = {\ frac {\ operatorname {Cov} _ { N-1} (X, Y) \ cdot (N-1) + {\ frac {N-1} {N}} (x_ {n} - {\ bar {x}} _ {n-1}) ( y_ {n} - {\ bar {y}} _ {n-1})} {N}} \\ & = {\ frac {\ operatorname {Cov} _ {N-1} (X, Y) \ cdot (N-1) + {\ frac {N} {N-1}} (x_ {n} - {\ bar {x}} _ {n}) (y_ {n} - {\ bar {y}} _ {n})} {N}}. \ end {alinhado}}}$

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)

Uma pequena modificação também pode ser feita 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)

Da mesma forma, existe uma fórmula para combinar as covariâncias de dois conjuntos que podem ser usados ​​para paralelizar o cálculo:

${\ displaystyle C_ {X} = C_ {A} + C_ {B} + ({\ bar {x}} _ {A} - {\ bar {x}} _ {B}) ({\ bar {y} } _ {A} - {\ bar {y}} _ {B}) \ cdot {\ frac {n_ {A} n_ {B}} {n_ {X}}}.}$

Versão ponderada em lote

Também existe uma versão do algoritmo online ponderado que atualiza em lote: vamos denotar os pesos e escrever ${\ displaystyle w_ {1}, \ dots w_ {N}}$

${\ displaystyle {\ begin {alignat} {2} {\ bar {x}} _ {n + k} & = {\ bar {x}} _ {n} & \, + \, & {\ frac {\ soma _ {i = n + 1} ^ {n + k} w_ {i} (x_ {i} - {\ bar {x}} _ {n})} {\ sum _ {i = 1} ^ {n + k} w_ {i}}} \\ {\ bar {y}} _ {n + k} & = {\ bar {y}} _ {n} & \, + \, & {\ frac {\ sum _ {i = n + 1} ^ {n + k} w_ {i} (y_ {i} - {\ bar {y}} _ {n})} {\ sum _ {i = 1} ^ {n + k} w_ {i}}} \\ C_ {n + k} & = C_ {n} & \, + \, & \ sum _ {i = n + 1} ^ {n + k} w_ {i} ( x_ {i} - {\ bar {x}} _ {n + k}) (y_ {i} - {\ bar {y}} _ {n}) \\ & = C_ {n} & \, + \ , & \ sum _ {i = n + 1} ^ {n + k} w_ {i} (x_ {i} - {\ bar {x}} _ {n}) (y_ {i} - {\ bar { y}} _ {n + k}) \\\ end {alinhado}}}$

A covariância pode então ser calculada como

${\ displaystyle \ operatorname {Cov} _ {N} (X, Y) = {\ frac {C_ {N}} {\ sum _ {i = 1} ^ {N} w_ {i}}}}$

Veja também

• Algoritmo de soma Kahan
• Desvios quadrados da média
• Método Yamartino

Referências

links externos

• Weisstein, Eric W. "Sample Variance Computation" . MathWorld .