Los algoritmos para calcular la varianza juegan un papel importante en la estadística computacional . Una dificultad clave en el diseño de buenos algoritmos para este problema es que las fórmulas para la varianza pueden involucrar sumas de cuadrados, lo que puede conducir a inestabilidad numérica así como a desbordamiento aritmético cuando se trabaja con 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 insesgada de la varianza de la población a partir de una muestra finita de n observaciones, la fórmula es:
Por lo tanto, un algoritmo ingenuo para calcular la varianza estimada viene dado por lo siguiente:
- Sea n ← 0, Sum ← 0, SumSq ← 0
- Para cada dato x :
- n ← n + 1
- Suma ← Suma + x
- SumSq ← SumSq + x × x
- Var = (SumSq - (Suma × Suma) / n) / (n - 1)
Este algoritmo se puede adaptar fácilmente para calcular la varianza de una población finita: simplemente divida por N en lugar de n - 1 en la última línea.
Debido a que SumSq y (Sum × Sum) / n pueden ser números muy similares, la cancelación puede llevar a que la precisión del resultado sea mucho menor que la precisión inherente de la aritmética de punto flotante utilizada para realizar el cálculo. Por tanto, este algoritmo no debería utilizarse en la práctica, [1] [2] y se han propuesto varios algoritmos alternativos numéricamente estables. [3] Esto es particularmente malo si la desviación estándar es pequeña en relación con la media. Sin embargo, el algoritmo se puede mejorar adoptando el método de la media supuesta .
Calcular datos desplazados
La varianza es invariante con respecto a los cambios en un parámetro de ubicación , una propiedad que se puede utilizar para evitar la cancelación catastrófica en esta fórmula.
con cualquier constante, que conduce a la nueva fórmula
cuanto más cerca es el valor medio, más preciso será el resultado, pero simplemente elegir un valor dentro del rango de muestras garantizará la estabilidad deseada. Si los valoresson pequeñas, 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 de la fórmula es siempre menor que el primero, por lo que no puede producirse ninguna cancelación. [2]
Si solo se toma la primera muestra como el algoritmo se puede escribir en el lenguaje de programación Python como
def shifted_data_variance ( data ): si len ( data ) < 2 : return 0.0 K = data [ 0 ] n = Ex = Ex2 = 0.0 para x en datos : n = n + 1 Ex + = x - K Ex2 + = ( x - K ) * ( x - K ) varianza = ( Ex2 - ( Ex * Ex ) / n ) / ( n - 1 ) # use n en lugar de (n-1) si desea calcular la varianza exacta de los datos dados # use (n-1) si los datos son muestras de una varianza de retorno de población mayor
Esta fórmula también facilita el cálculo incremental que se puede expresar como
K = n = Ex = Ex2 = 0.0def add_variable ( x ): global K , n , Ex , Ex2 si 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 / ndef get_variance (): n global , Ex , Ex2 return ( Ex2 - ( Ex * Ex ) / n ) / ( n - 1 )
Algoritmo de dos pasos
Un enfoque alternativo, utilizando una fórmula diferente para la varianza, primero calcula la media de la muestra,
y luego calcula la suma de los cuadrados de las diferencias de la media,
donde s es la desviación estándar. Esto viene dado por el siguiente código:
def two_pass_variance ( datos ): n = sum1 = sum2 = 0 para x en los datos : n + = 1 suma1 + = x media = suma1 / n para x en datos : suma2 + = ( x - media ) * ( x - media ) varianza = suma2 / ( n - 1 ) devolver varianza
Este algoritmo es numéricamente estable si n es pequeño. [1] [4] Sin embargo, los resultados de estos dos algoritmos simples ("ingenuo" y "de dos pasadas") pueden depender excesivamente del orden de los datos y pueden dar resultados deficientes para conjuntos de datos muy grandes debido al redondeo repetido Error en la acumulación de las sumas. Se pueden utilizar técnicas como la suma compensada para combatir este error hasta cierto punto.
Algoritmo en línea de Welford
A menudo es útil poder calcular la varianza en una sola pasada , inspeccionando cada valorsólo una vez; por ejemplo, cuando los datos se recopilan sin suficiente almacenamiento para mantener todos los valores, o cuando los costos de acceso a la memoria dominan los de la computación. Para tal algoritmo en línea , se requiere una relación de recurrencia entre cantidades a partir de las cuales se pueden calcular las estadísticas requeridas de una manera 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 x n . Aquí, x n denota la media muestral de las primeras n muestras ( x 1 , ..., x n ), s2
nsu varianza muestral, y σ2
n 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 escala con n . Una mejor cantidad para actualizar es la suma de cuadrados de las diferencias de la media actual,, aquí denotado :
Este algoritmo fue encontrado por Welford, [5] [6] y ha sido analizado a fondo. [2] [7] También es común denotar y . [8]
A continuación se proporciona un ejemplo de implementación de Python para el algoritmo de Welford.
# Para un nuevo valor newValue, calcule el nuevo recuento, la nueva media, el nuevo M2. # mean acumula la media de todo el conjunto de datos # M2 agrega la distancia al cuadrado de la media # count agrega el número de muestras vistas hasta ahora def update ( existingAggregate , newValue ): ( count , mean , M2 ) = existingAggregate count + = 1 delta = newValue - media media + = delta / count delta2 = newValue - mean M2 + = delta * delta2 return ( count , mean , M2 )# Recuperar la media, la varianza y la varianza de la muestra a partir de un agregado def finalize ( existingAggregate ): ( recuento , media , M2 ) = existingAggregate si el recuento < 2 : el retorno del flotador ( "nan" ) otra cosa : ( media , varianza , sampleVariance ) = ( mean , M2 / count , M2 / ( count - 1 )) return ( media , varianza , muestraVariance )
Este algoritmo es mucho menos propenso a perder precisión debido a una cancelación catastrófica , pero podría no ser tan eficiente debido a la operación de división dentro del bucle. Para un algoritmo de dos pasos particularmente robusto para calcular la varianza, primero se puede calcular y restar una estimación de la media y luego usar este algoritmo en los residuos.
El algoritmo paralelo a continuación ilustra cómo combinar varios conjuntos de estadísticas calculadas en línea.
Algoritmo incremental ponderado
El algoritmo se puede ampliar para manejar pesos de muestra desiguales, reemplazando el contador simple n con la suma de pesos vistos hasta ahora. West (1979) [9] sugiere este algoritmo incremental :
def weighted_incremental_variance ( data_weight_pairs ): w_sum = w_sum2 = media = S = 0 para x , w en data_weight_pairs : # Alternativamente "para x, w en zip (datos, pesos):" 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 ) población_varianza = S / w_sum # Corrección de Bessel para muestras ponderadas # Ponderaciones de frecuencia sample_frequency_variance = S / ( w_sum - 1 ) # Ponderaciones de confiabilidad sample_reliability_variance = S / ( w_sum - w_sum2 / w_sum )
Algoritmo paralelo
Chan y col. [10] tenga en cuenta que el algoritmo en línea de Welford detallado anteriormente es un caso especial de un algoritmo que funciona para combinar conjuntos arbitrarios y :
- .
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 lo es en el caso. En tales casos, prefiera.
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
Esto se puede generalizar para permitir la paralelización con AVX , con GPU y clústeres de computadoras , y con covarianza. [3]
Ejemplo
Suponga que todas las operaciones de coma flotante utilizan aritmética de doble precisión estándar IEEE 754 . Considere la muestra (4, 7, 13, 16) de una población infinita. Con base en esta muestra, la media poblacional estimada es 10 y la estimación insesgada de la varianza poblacional es 30. Tanto el algoritmo ingenuo como el algoritmo de dos pasos calculan estos valores correctamente.
A continuación, considere la muestra ( 10 8 + 4 , 10 8 + 7 , 10 8 + 13 , 10 8 + 16 ), que da lugar a la misma varianza estimada que la primera muestra. El algoritmo de dos pasos calcula esta estimación de la varianza correctamente, pero el algoritmo ingenuo devuelve 29.333333333333332 en lugar de 30.
Si bien esta pérdida de precisión puede ser tolerable y verse como un defecto menor del algoritmo ingenuo, aumentar aún más el desplazamiento hace que el error sea catastrófico. Considere la muestra ( 10 9 + 4 , 10 9 + 7 , 10 9 + 13 , 10 9 + 16 ). Una vez más, la varianza de la población estimada de 30 se calcula correctamente mediante el algoritmo de dos pasos, pero el algoritmo ingenuo ahora la calcula como −170,66666666666666. Este es un problema grave con el algoritmo ingenuo y se debe a una cancelación catastrófica en la resta de dos números similares en la etapa final del algoritmo.
Estadísticas de orden superior
Terriberry [11] extiende las fórmulas de Chan para calcular el tercer y cuarto momento central , necesario, por ejemplo, al estimar la asimetría y la curtosis :
Aquí el son de nuevo las sumas de potencias de las diferencias de la media , donación
Para el caso incremental (es decir, ), esto se simplifica a:
Preservando 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 costo incremental mínimo.
Un ejemplo del algoritmo en línea para la curtosis implementado como se describe es:
def kurtosis_en línea ( datos ): n = media = M2 = M3 = M4 = 0 para x en datos : n1 = n n = n + 1 delta = x - media delta_n = delta / n delta_n2 = delta_n * delta_n term1 = delta * delta_n * n1 media = media + delta_n M4 = M4 + term1 * delta_n2 * ( n * n - 3 * n + 3 ) + 6 * delta_n2 * M2 - 4 * delta_n * M3 M3 = M3 + término1 * delta_n * ( n - 2 ) - 3 * delta_n * M2 M2 = M2 + término1 # Tenga en cuenta que también puede calcular la varianza usando M2 y la asimetría usando M3 # Precaución: Si todas las entradas son iguales, M2 será 0, lo que dará como resultado una división entre 0. curtosis = ( n * M4 ) / ( M2 * M2 ) - 3 curtosis de retorno
Pébaÿ [12] extiende estos resultados a momentos centrales de orden arbitrario , para los casos incrementales y por pares, y posteriormente Pébaÿ et al. [13] para momentos ponderados y compuestos. También se pueden encontrar fórmulas similares para la covarianza .
Choi y Sweetman [14] ofrecen dos métodos alternativos para calcular la asimetría y la curtosis, cada uno de los cuales puede ahorrar importantes requisitos de memoria de la computadora y tiempo de CPU en ciertas aplicaciones. El primer enfoque consiste en calcular los momentos estadísticos separando los datos en contenedores y luego calculando los momentos a partir de la geometría del histograma resultante, que efectivamente se convierte en un algoritmo de un solo paso para momentos superiores. Un beneficio es que los cálculos estadísticos de momento se pueden llevar a cabo con precisión arbitraria, de modo que los cálculos se pueden sintonizar con la precisión de, por ejemplo, el formato de almacenamiento de datos o el hardware de medición original. Se puede construir un histograma relativo de una variable aleatoria de la manera convencional: el rango de valores potenciales se divide en contenedores y el número de ocurrencias dentro de cada contenedor se cuenta y se traza de manera que el área de cada rectángulo sea igual a la porción de los valores de la muestra. dentro de ese contenedor:
dónde y representar la frecuencia y la frecuencia relativa en bin y es el área total del histograma. Después de esta normalización, el momentos crudos y momentos centrales de se puede calcular a partir del histograma relativo:
donde el superíndice indica que los momentos se calculan a partir del histograma. Para un ancho de contenedor constante estas dos expresiones se pueden simplificar usando :
El segundo enfoque de Choi y Sweetman [14] es una metodología analítica para combinar momentos estadísticos de segmentos individuales de una historia de tiempo de modo que los momentos generales resultantes sean los de la historia de tiempo completa. Esta metodología podría utilizarse para el cálculo paralelo de momentos estadísticos con la combinación posterior de esos momentos, o para la combinación de momentos estadísticos calculados en tiempos secuenciales.
Si Se conocen conjuntos de momentos estadísticos: por , luego cada se puede expresar en términos del equivalente momentos crudos:
dónde se toma generalmente como la duración de la historial de tiempo, o el número de puntos si es constante.
El beneficio de expresar los momentos estadísticos en términos de es que el Los conjuntos se pueden combinar por suma, y no hay límite superior en el valor de .
donde el subíndice representa la historia de tiempo concatenada o combinada . Estos valores combinados de luego se puede transformar inversamente en momentos crudos que representan la historia del tiempo concatenada completa
Relaciones conocidas entre los momentos crudos () y los momentos centrales () se utilizan para calcular los momentos centrales de la historia de tiempo concatenada. Finalmente, los momentos estadísticos de la historia concatenada se calculan a partir de los momentos centrales:
Covarianza
Se pueden utilizar algoritmos muy similares para calcular la covarianza .
Algoritmo ingenuo
El algoritmo ingenuo es:
Para el algoritmo anterior, se podría usar el siguiente código de Python:
def naive_covariance ( data1 , data2 ): n = len ( data1 ) sum12 = 0 sum1 = sum ( data1 ) sum2 = sum ( data2 ) para i1 , i2 en zip ( data1 , data2 ): sum12 + = i1 * i2 covarianza = ( suma12 - suma1 * suma2 / n ) / n devuelve covarianza
Con estimación de la media
En cuanto a la varianza, la covarianza de dos variables aleatorias también es invariante al desplazamiento, por lo que dados dos valores constantes y se puede escribir:
y, de nuevo, la elección de un valor dentro del rango de valores estabilizará la fórmula frente a cancelaciones catastróficas y la hará más robusta frente a grandes sumas. Tomando el primer valor de cada conjunto de datos, el algoritmo se puede escribir como:
def covarianza_datos_desplazados ( datos_x , datos_y ): n = len ( datos_x ) si n < 2 : return 0 kx = datos_x [ 0 ] ky = datos_y [ 0 ] Ex = Ey = Exy = 0 para ix , iy en zip ( datos_x , datos_y ): Ex + = ix - kx Ey + = iy - ky Exy + = ( ix - kx ) * ( iy - ky ) return ( Exy - Ex * Ey / n ) / n
Dos pasadas
El algoritmo de dos pasos primero calcula las medias de la muestra y luego la covarianza:
El algoritmo de dos pasos se puede escribir como:
def two_pass_covariance ( dato1 , dato2 ): n = len ( dato1 ) media1 = suma ( datos1 ) / n media2 = suma ( datos2 ) / n covarianza = 0 para i1 , i2 en zip ( data1 , data2 ): a = i1 - mean1 b = i2 - mean2 covariance + = a * b / n return covariance
Una versión compensada un poco más precisa realiza el algoritmo ingenuo completo en los residuos. Las sumas finales y debe ser cero, pero la segunda pasada compensa cualquier pequeño error.
En línea
Existe un algoritmo estable de un solo paso, similar al algoritmo en línea para calcular la varianza, que calcula el co-momento :
La aparente asimetría en esa última ecuación se debe al hecho de que , por lo que ambos términos de actualización son iguales a . Se puede lograr una precisión aún mayor calculando primero los medios y luego utilizando el algoritmo estable de una pasada en los residuos.
Por tanto, la covarianza se puede calcular como
def online_covariance ( datos1 , datos2 ): MEANX = meany = C = n = 0 para x , y en zip ( data1 , data2 ): n + = 1 dx = x - MEANX MEANX + = dx / n meany + = ( y - media ) / n C + = dx * ( y - media ) población_covar = C / n # Corrección de Bessel para la varianza de la muestra sample_covar = C / ( n - 1 )
También se puede realizar una pequeña modificación para calcular la covarianza ponderada:
def online_weighted_covariance ( datos1 , datos2 , Data3 ): MEANX = meany = 0 wsum = wsum2 = 0 C = 0 para x , y , w en zip ( datos1 , datos2 , Data3 ): wsum + = w wsum2 + = w * w dx = x - mediax mediax + = ( w / wsum ) * dx meany + = ( w / wsum ) * ( y - meany ) C + = w * dx * ( y - meany ) población_covar = C / wsum # Corrección de Bessel para la varianza de la muestra # Ponderaciones de frecuencia sample_frequency_covar = C / ( wsum - 1 ) # Ponderaciones de confiabilidad sample_reliability_covar = C / ( wsum - wsum2 / wsum )
Asimismo, existe una fórmula para combinar las covarianzas de dos conjuntos que se puede utilizar para paralelizar el cálculo: [3]
Versión por lotes ponderada
También existe una versión del algoritmo en línea ponderado que se actualiza por lotes: let denotar los pesos y escribir
La covarianza se puede calcular como
Ver también
- Algoritmo de suma de Kahan
- Desviaciones cuadradas de la media
- Método Yamartino
Referencias
- ↑ a b Einarsson, Bo (2005). Precisión y confiabilidad en la informática científica . SIAM. pag. 47. ISBN 978-0-89871-584-2.
- ^ a b c Chan, Tony F .; Golub, Gene H .; LeVeque, Randall J. (1983). "Algoritmos para calcular la varianza muestral: análisis y recomendaciones" (PDF) . El estadístico estadounidense . 37 (3): 242–247. doi : 10.1080 / 00031305.1983.10483115 . JSTOR 2683386 .
- ^ a b c Schubert, Erich; Gertz, Michael (9 de julio de 2018). Cálculo paralelo numéricamente estable de (co) varianza . ACM. pag. 10. doi : 10.1145 / 3221269.3223036 . ISBN 9781450365055. S2CID 49665540 .
- ^ Higham, Nicholas (2002). Precisión y estabilidad de los algoritmos numéricos (2 ed) (Problema 1.10) . SIAM.
- ^ Welford, BP (1962). "Nota sobre un método para calcular sumas corregidas de cuadrados y productos". Tecnometría . 4 (3): 419–420. doi : 10.2307 / 1266577 . JSTOR 1266577 .
- ^ Donald E. Knuth (1998). El arte de la programación informática , volumen 2: Algoritmos seminuméricos , 3ª ed., P. 232. Boston: Addison-Wesley.
- ^ Ling, Robert F. (1974). "Comparación de varios algoritmos para calcular medias y varianzas muestrales". Revista de la Asociación Estadounidense de Estadística . 69 (348): 859–866. doi : 10.2307 / 2286154 . JSTOR 2286154 .
- ^ http://www.johndcook.com/standard_deviation.html
- ^ West, DHD (1979). "Actualización de las estimaciones de la media y la varianza: un método mejorado". Comunicaciones de la ACM . 22 (9): 532–535. doi : 10.1145 / 359146.359153 . S2CID 30671293 .
- ^ Chan, Tony F .; Golub, Gene H .; LeVeque, Randall J. (1979), "Actualización de fórmulas y un algoritmo por pares para calcular varianzas muestrales". (PDF) , Informe técnico STAN-CS-79-773 , Departamento de Ciencias de la Computación, Universidad de Stanford.
- ^ Terriberry, Timothy B. (2007), Computing Higher-Order Moments Online , archivado desde el original el 23 de abril de 2014 , consultado el 5 de mayo de 2008
- ^ Pébaÿ, Philippe (2008), "Fórmulas para el cálculo paralelo robusto de covarianzas y momentos estadísticos de orden arbitrario" (PDF) , Informe técnico SAND2008-6212 , Sandia National Laboratories[ enlace muerto permanente ]
- ^ Pébaÿ, Philippe; Terriberry, Timothy; Kolla, Hemanth; Bennett, Janine (2016), "Fórmulas escalables y numéricamente estables para el cálculo paralelo y en línea de momentos centrales multivariados de orden superior con pesos arbitrarios" , Estadística computacional , Springer, 31 (4): 1305-1325, doi : 10.1007 / s00180- 015-0637-z , S2CID 124570169
- ^ a b Choi, Myoungkeun; Sweetman, Bert (2010), "Efficient Calculation of Statistical Moments for Structural Health Monitoring", Journal of Structural Health Monitoring , 9 (1): 13-24, doi : 10.1177 / 1475921709341014 , S2CID 17534100
enlaces externos
- Weisstein, Eric W. "Ejemplo de cálculo de la varianza" . MathWorld .