Algoritmos para cálculo da variância desempenham um papel importante na estatística computacional. Uma dificuldade central no projeto de bons algoritmos para esse problema é que as fórmulas da variância podem envolver somas de quadrados, o que pode levar à instabilidade numérica, bem como a 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 é:
Usando a correção de Bessel para calcular uma estimativa não enviesada da variância populacional a partir de uma amostra finita de n observações, a fórmula é:
Portanto, um algoritmo ingênuo para calcular a variância estimada é dado por:
Seja Para cada dado :
Esse 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.
Como e 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. Assim, esse algoritmo não deve ser usado na prática, e foram propostos diversos algoritmos alternativos numericamente estáveis. Esse problema é particularmente grave quando o desvio padrão é pequeno em relação à média.
Cálculo com dados deslocados
A variância é invariante em relação a alterações de um parâmetro de localização, propriedade que pode ser usada para evitar o cancelamento catastrófico nessa fórmula.
com sendo uma constante qualquer, o que leva à nova fórmula
Quanto mais próximo estiver do valor médio, mais preciso será o resultado, mas escolher qualquer valor dentro do intervalo da amostra já fornece a estabilidade desejada. Se os valores forem pequenos, não haverá problemas com a soma de seus quadrados; por outro lado, se forem grandes, isso necessariamente significa que a variância também é grande. Em qualquer caso, o segundo termo da fórmula é sempre menor que o primeiro; portanto, não pode ocorrer cancelamento.
O algoritmo pode ser escrito na linguagem de programação Python como def shifted_data_variance(data, K): if len(data) < 2: return 0.0 n = Ex = Ex2 = 0.0 for x in data: n += 1 Ex += x - K Ex2 += (x - K) ** 2 variance = (Ex2 - Ex**2 / n) / (n - 1) # use n em vez de (n-1) para calcular a variância exata dos dados fornecidos # use (n-1) se os dados forem amostras de uma população maior return variance
Algoritmo de duas passagens Uma abordagem alternativa, que usa uma fórmula diferente para a variância, primeiro calcula a média amostral,
e então calcula a soma dos quadrados das diferenças em relação à média,
em que s é o desvio padrão. Isso é dado pelo seguinte código:
def two_pass_variance(data): n = len(data) mean = sum(data) / n variance = sum((x - mean) ** 2 for x in data) / (n - 1) return variance
Esse trecho de código, quando executado no CPython 3.12 ou posterior, é sempre numericamente estável. Isso ocorre porque essas versões usam o esquema de soma compensada de Neumaier na função sum(), tornando-a resistente a erros repetidos de arredondamento. Muitas linguagens voltadas à computação numérica oferecem construções semelhantes.
Esse algoritmo, se implementado com uma soma ingênua sum = 0; for x in data: sum += x; mean = sum / n, seria numericamente estável somente se n fosse pequeno, por causa de erros de arredondamento repetidos ao adicionar um número pequeno a outro muito, muito maior. Metadados também listados em ACM Digital Library.
Algoritmos incrementais Frequentemente é útil poder calcular a variância em uma única passagem, examinando 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 predominam sobre os custos de computação. Para esse tipo de algoritmo online, é necessária uma relação de recorrência entre grandezas a partir das quais as estatísticas desejadas possam ser calculadas de modo numericamente estável.
Dados deslocados incrementais O algoritmo de dados deslocados apresentado anteriormente inclui a atualização simples dos dados em um laço, o que o torna naturalmente incremental. Ele pode ser expresso como:
class ShiftDataVariance: def __init__(self): self.K = 0.0 self.n = 0 self.Ex = 0.0 self.Ex2 = 0.0
def add_variable(self, x: float): if self.n == 0: self.K = x self.n += 1 self.Ex += x - self.K self.Ex2 += (x - self.K) ** 2
def remove_variable(self, x: float): self.n -= 1 self.Ex -= x - self.K self.Ex2 -= (x - self.K) ** 2
def add_variables(self, xs: list[float]): # Python usa o algoritmo de Neumaier na função embutida sum(), # que é mais preciso do que um laço simples. if self.n == 0 and xs: self.K = xs self.n += len(xs) self.Ex += sum(x - self.K for x in xs) self.Ex2 += sum((x - self.K) ** 2 for x in xs)
def get_mean(self) -> float: return self.K + self.Ex / self.n
def get_variance(self) -> float: return (self.Ex2 - self.Ex**2 / self.n) / (self.n - 1)
Algoritmo online de Welford De forma análoga ao algoritmo de duas passagens acima, continua sendo desejável usar a média real dos dados, em vez de simplesmente escolher o primeiro elemento. As fórmulas a seguir podem ser usadas para atualizar a média e a variância (estimada) da sequência com um elemento adicional xn. Aqui, denota a média amostral das primeiras n amostras , sua variância amostral enviesada, e sua variância amostral não enviesada.
Essas fórmulas sofrem de instabilidade numérica, pois subtraem repetidamente um número pequeno de um número grande que cresce com n. Uma quantidade melhor para atualização é a soma dos quadrados das diferenças em relação à média atual, , aqui denotada por :
O algoritmo a seguir foi desenvolvido por Welford,Donald E. Knuth (1998). The Art of Computer Programming, volume 2: Seminumerical Algorithms, 3.ª ed., p. 232. Boston: Addison-Wesley. e foi analisado em profundidade. Também é comum denotar e . A seguir é apresentada uma implementação de exemplo em Python para o algoritmo de Welford, usando a mesma estrutura do algoritmo de "dados deslocados" acima:
class WelfordVariance: def __init__(self): # Comparação com ShiftDataVariance: self.mean = 0.0 # = K + Ex / n self.count = 0 # = n self.M2 = 0.0 # = Ex2 - (Ex)^2 / n
def add_variable(self, x: float): self.count += 1 old_mean = self.mean self.mean += (x - self.mean) / self.count self.M2 += (x - old_mean) * (x - self.mean)
def remove_variable(self, x: float): self.count -= 1 new_mean = self.mean self.mean -= (x - self.mean) / self.count self.M2 -= (x - new_mean) * (x - self.mean)
def get_mean(self) -> float: return self.mean
def get_variance(self) -> float: return self.M2 / self.count
def get_sample_variance(self) -> float: return self.M2 / (self.count - 1)
Esse algoritmo é muito menos propenso à perda de precisão decorrente de cancelamento catastrófico, mas pode não ser tão eficiente devido à operação de divisão dentro do laço. Para um algoritmo de duas passagens particularmente robusto no cálculo da variância, pode-se primeiro calcular e subtrair uma estimativa da média e, em seguida, aplicar esse algoritmo aos resíduos.
O algoritmo paralelo abaixo ilustra como combinar múltiplos 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 observados até então. West (1979) propõe este algoritmo incremental:
from collections import namedtuple WeightedVariances = namedtuple("WeightedVariances", "pop freq reli") 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**2 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 # Correção de Bessel para amostras ponderadas # Pesos de frequência sample_frequency_variance = S / (w_sum - 1)
# Pesos de confiabilidade sample_reliability_variance = S / (w_sum - w_sum2 / w_sum) return WeightedVariances(population_variance, sample_frequency_variance, sample_reliability_variance)
Algoritmo paralelo Chan et al. observam que o algoritmo online de Welford detalhado acima é um caso particular de um algoritmo que permite combinar conjuntos arbitrários e : . Isso pode ser útil quando, por exemplo, várias unidades de processamento são atribuídas a partes distintas 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 mesma forma que no caso . Nesses casos, prefira . Reutilizando a classe WelfordVariance acima, temos:
def merge(a: WelfordVariance, b: WelfordVariance) -> WelfordVariance: ab = WelfordVariance() ab.count = a.count + b.count delta = b.mean - a.mean ab.mean = (a.count * a.mean + b.count * b.mean) / ab.count ab.M2 = a.M2 + b.M2 + delta**2 * a.count * b.count / ab.count return ab
# exemplo: ab = merge(a, b) print(ab.get_sample_variance())
Esse algoritmo permite dividir um conjunto de dados em várias partes, processá-las em paralelo e, em seguida, combinar os resultados. Isso possibilita paralelização de várias formas, inclusive com AVX, GPUs e clusters de computadores. O algoritmo também pode ser adaptado para estatísticas de ordem superior e para covariância.
Exemplo Suponha que todas as operações de ponto flutuante usem a aritmética padrão IEEE 754 de precisão dupla. Considere a amostra (4, 7, 13, 16) de uma população infinita. Com base nessa amostra, a média populacional estimada é 10 e a estimativa não enviesada da variância populacional é 30. Tanto o algoritmo ingênuo quanto o algoritmo de duas passagens calculam corretamente esses valores.
Considere agora a amostra (, , , ), que produz a mesma variância estimada da primeira amostra. O algoritmo de duas passagens calcula corretamente essa estimativa de variância, 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 pequena falha do algoritmo ingênuo, aumentar ainda mais o deslocamento torna o erro catastrófico. Considere a amostra (, , , ). Novamente, a variância populacional estimada de 30 é calculada corretamente pelo algoritmo de duas passagens, mas o algoritmo ingênuo agora a calcula como −170.66666666666666. Esse é um problema grave do algoritmo ingênuo e decorre do cancelamento catastrófico na subtração de dois números semelhantes na etapa final do algoritmo.
Estatísticas de ordem superior Terriberry estende as fórmulas de Chan ao cálculo do terceiro e do quarto momentos centrais, necessários, por exemplo, para estimar assimetria e curtose:
Aqui, os são novamente as somas das potências das diferenças em relação à média, , resultando em
Para o caso incremental (isto é, ), isso se simplifica para:
Ao preservar o valor , é necessária apenas uma operação de divisão e, assim, as estatísticas de ordem superior podem ser calculadas com baixo custo incremental.
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**2 term1 = delta * delta_n * n1 mean = mean + delta_n M4 = M4 + term1 * delta_n2 * (n**2 - 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
# Observação: também é possível calcular a variância usando M2 e a assimetria usando M3 # Atenção: se todas as entradas forem iguais, M2 será 0, resultando em divisão por 0. kurtosis = (n * M4) / (M2**2) - 3 return kurtosis
Pébaÿ estende ainda mais esses resultados para momentos centrais de ordem arbitrária, nos casos incremental e par a par, e posteriormente Pébaÿ et al.