Los algoritmos para calcular la varianza juegan un papel importante en la estadística informática. Una dificultad clave en el diseño de buenos algoritmos para este problema es que las fórmulas de la varianza pueden incluir sumas de cuadrados, lo que puede llevar a problemas de inestabilidad numérica así como a rebosamiento cuando se trata de valores grandes.
Algoritmo "ingenuo" Una fórmula para calcular la varianza de una población completa de tamaño N es:
Usando la corrección de Bessel para calcular una estimación no sesgada de la varianza poblacional de una muestra finita de observaciones n, la fórmula es:
Por lo tanto, un algoritmo "ingenuo" para calcular la varianza estimada viene dado por el procedimiento:
Let For each datum :
Este algoritmo se puede adaptar fácilmente para calcular la varianza de una población finita: simplemente, basta con dividr entre N en lugar de entre n - 1 en la última línea.
Como y pueden ser números muy similares, cancelación de los restos puede llevar a que la precisión del resultado sea mucho menor que la precisión inherente de coma flotante utilizada para realizar el cálculo. Por lo tanto, se ha propuesto este algoritmo no debe utilizarse en la práctica. y varios algoritmos alternativos, numéricamente estables. Esto es particularmente malo si la desviación estándar es pequeña en relación con la media. Sin embargo, el algoritmo puede mejorarse adoptando el método de la media falsa.
Cálculo de datos desplazados
Podemos usar una propiedad de la varianza para evitar la cancelación catastrófica en esta fórmula, a saber, la varianza es invariante con respecto a los cambios respecto a un parámetro de localización
Siendo una constante cualquiera, se obtiene la nueva fórmula:
Cuanto más cerca esté del valor medio, más preciso será el resultado, pero simplemente eligiendo un valor dentro del rango de muestras se garantizará la estabilidad deseada. Si los valores son pequeños, entonces no hay problemas con la suma de sus cuadrados, por el contrario, si son grandes, necesariamente significa que la varianza también es grande. En cualquier caso, el segundo término en la fórmula es siempre más pequeño que el primero, por lo que no puede producirse ninguna cancelación.
Si se toma solo la primera muestra como , el algoritmo se puede escribir en el lenguaje de programación Python como
def shifted_data_variance(data): if len(data) < 2: return 0.0 K = data 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 también el cálculo incremental, que puede expresarse como K = n = Ex = Ex2 = 0.0
def add_variable (x): si (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 (): volver K + Ex / n
def get_variance (): retorno (Ex2 - (Ex * Ex) / n) / (n-1)
Algoritmo de dos pasos Una aproximación alternativa, usando una fórmula diferente para la varianza, computa primero la media de la muestra,
y después computa la suma de los cuadrados de las diferencias con la media,
donde s es la desviación estándar. Este planteamiento viene dado por el siguiente 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 es numéricamente estable si n es pequeño. Sin embargo, los resultados de estos dos algoritmos simples ("ingenuo" y "de dos pasos") pueden depender excesivamente del orden de los datos y pueden dar resultados deficientes para conjuntos de datos muy grandes, debido a repetidos errores de redondeo en la acumulación de las sumas. Técnicas como la suma compensada se pueden usar para combatir este error hasta cierto punto.
Algoritmo en línea de Welford A menudo es útil poder calcular la variación en una sola pasada, inspeccionando cada valor solo una vez; por ejemplo, cuando los datos se recopilan sin el almacenamiento suficiente para mantener todos los valores, o cuando los costos de acceso a la memoria dominan a los de cálculo. Para este algoritmo en línea, se utiliza una relación de recurrencia entre cantidades a partir de las cuales las estadísticas requeridas se pueden calcular de forma numéricamente estable.
Las siguientes fórmulas se pueden utilizar para actualizar la media y la varianza (estimada) de la secuencia, para un elemento adicional xn. Aquí, n denota la media muestral de las primeras n muestras (x1, ..., xn), s su varianza muestral, y σ su varianza poblacional.
Estas fórmulas sufren de inestabilidad numérica, ya que restan repetidamente un número pequeño de un número grande que se escala con n. Una mejor cantidad para actualizar es la suma de cuadrados de las diferencias de la media actual, , que aquí denota :
Este algoritmo fue ideado por Welford,Donald Knuth (1998). The Art of Computer Programming, volume 2: Seminumerical Algorithms, 3rd edn., p. 232. Boston: Addison-Wesley. y se ha analizado exhaustivamente. También es común denotar y .https://www.johndcook.com/standard_deviation.html
A continuación se muestra un ejemplo de implementación en Python para el algoritmo de Welford.
# 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 es mucho menos propenso a la pérdida de precisión debido a la cancelación catastrófica, pero podría no ser tan eficiente debido a la operación de división incluida dentro del bucle. Para disponer de un algoritmo de dos pasos particularmente robusto con el fin de calcular la varianza, primero se puede calcular y restar una estimación de la media y luego usar este algoritmo con los residuos.
El algoritmo paralelo mostrado más adelante ilustra cómo combinar varios conjuntos de estadísticas calculadas en línea.
Algoritmo incremental ponderado El algoritmo puede extenderse para manejar pesos de muestra desiguales, reemplazando el simple contador n con la suma de pesos vista hasta ahora. West (1979) sugiere este algoritmo incremental:
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.. observaron que el algoritmo en línea de Welford detallado anteriormente es un caso especial de un algoritmo que funciona para cualquier partición de la muestra en los conjuntos , :
. Esto puede ser útil cuando, por ejemplo, se pueden asignar múltiples unidades de procesamiento a partes discretas de la entrada.
El método de Chan para estimar la media es numéricamente inestable cuando y ambos son grandes, porque el error numérico en no se reduce de la forma en que se encuentra en el caso . En tales casos, es preferible .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) Esto se puede generalizar para permitir la paralelización con Extensiones Vectoriales Avanzadas (EVA), con unidades de procesamiento gráfico y agregados de computadoras, y para la covarianza.
Ejemplo Supóngase que todas las operaciones de punto flotante utilizan aritmética IEEE 754 de doble precisión estándar. Considérese la muestra (4, 7, 13, 16) de una población infinita. Basándose en esta muestra, la media poblacional estimada es 10, y la estimación no sesgada de la varianza poblacional es 30. Tanto el algoritmo "ingenuo" como el algoritmo de dos pasos calculan estos valores correctamente.
A continuación, considérese la muestra (, , , ), que da lugar a la misma varianza estimada que la primera muestra. El algoritmo de dos pasos calcula esta estimación de varianza correctamente, pero el algoritmo "ingenuo" devuelve 29.33333333333333332 en lugar de 30.
Si bien esta pérdida de precisión puede ser tolerable y vista como un defecto menor del algoritmo "ingenuo", un aumento adicional del desplazamiento hace que el error sea catastrófico. Considérese la muestra (, , , ). De nuevo, la varianza estimada de la población de 30 se calcula correctamente mediante el algoritmo de dos pasos, pero el algoritmo "ingenuo" ahora lo calcula como −170.66666666666666. Este es un problema serio con el algoritmo "ingenuo" y se debe a la cancelación catastrófica en la resta de dos números similares en la etapa final del algoritmo.
Estadísticas de orden superior Terriberry amplía las fórmulas de Chan para calcular el tercer y cuarto momento central, necesarios, por ejemplo, al estimar la asimetría estadística y la curtosis:
Aquí, los son nuevamente las sumas de las potencias de las diferencias conrespecto a la media , dando
Para el caso incremental (es decir, ), esto se simplifica a:
Al preservar el valor , solo se necesita una operación de división y, por lo tanto, las estadísticas de orden superior se pueden calcular con un reducido costo incremental.
Un ejemplo del algoritmo en línea para la curtosis implementado como se describe es: 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ÿ extiende aún más estos resultados a momentos centrales de orden arbitrario, para los casos incrementales y por pares, y posteriormente Pébaÿ et al. lo hicieron para momentos ponderados y compuestos. También se pueden encontrar fórmulas similares para la covarianza.
Choi y Sweetman