Algoritmos para calcular la varianza

Compartir Imprimir Citar

algoritmos importantes en estadísticas numéricas

Los algoritmos para calcular la varianza juegan un papel importante en las estadísticas computacionales. Una dificultad clave en el diseño de buenos algoritmos para este problema es que las fórmulas para la varianza pueden incluir sumas de cuadrados, lo que puede provocar inestabilidad numérica y desbordamiento aritmético cuando se trata de valores grandes.

Algoritmo ingenuo

Una fórmula para calcular la varianza de una población entera de tamaño N es:

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

Usando la corrección de Bessel para calcular una estimación imparcial de la varianza de la población a partir de una muestra finita de n observaciones, la fórmula es:

{displaystyle s^{2}=left({frac {sum _{i=1}^{n}x_{i}^{2}}{n}}-left({frac {sum _{i=1}^{n}x_{i}}{n}}right)^{2}right)cdot {frac {n}{n-1}}.}

Por lo tanto, un algoritmo ingenuo para calcular la varianza estimada viene dado por lo siguiente:

Vamos $n \leftarrow 0, Sum \leftarrow 0, SumSq \leftarrow 0$
Para cada dato x:
- $n \leftarrow n + 1$
- $Sum \leftarrow Sum + x$
- $SumSq \leftarrow SumSq + x \times x$
$Var = (SumSq - (Sum \times Sum) / (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\timesSum)/ n$ pueden ser números muy similares, la cancelación puede hacer que la precisión del resultado sea mucho menor que la precisión inherente de la aritmética de coma flotante utilizada para realizar el cálculo. Por lo tanto, este algoritmo no debe usarse en la práctica, y se han propuesto 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.

Cómputo de datos desplazados

La varianza es invariable con respecto a los cambios en un parámetro de ubicación, una propiedad que se puede usar para evitar la cancelación catastrófica en esta fórmula.

{displaystyle operatorname {Var} (X-K)=operatorname {Var} (X).}

con $K$ cualquier constante, que conduce a la nueva fórmula

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

más cerca $K$ es al valor medio el resultado será más exacto, pero sólo elegir un valor dentro del rango de muestras garantizará la estabilidad deseada. Si los valores $(x_i - K)$ son pequeños entonces no hay problemas con la suma de sus cuadrados, por el contrario, si son grandes necesariamente significa que la varianza es grande también. En cualquier caso el segundo término en la fórmula es siempre menor que el primero por lo tanto no puede ocurrir cancelación.

Si sólo la primera muestra se toma como $K$ el algoritmo se puede escribir en Python lenguaje de programación como

def cambio_data_varianza()datos): si Len()datos) . 2: retorno 0,0 K = datos[0] n = Ex = Ex2 = 0,0 para x dentro datos: n += 1 Ex += x - K Ex2 += ()x - K)#2 diferencia = ()Ex2 - Ex#2 / n) / ()n - 1) # use n en lugar de (n-1) si desea calcular la varianza exacta de los datos dados # use (n-1) if data are samples of a larger population retorno diferencia

Esta fórmula también facilita el cálculo incremental que se puede expresar como

K = Ex = Ex2 = 0,0n = 0def add_variable()x): mundial K, n, Ex, Ex2 si n == 0: K = x n += 1 Ex += x - K Ex2 += ()x - K)#2def remove_variable()x): mundial K, n, Ex, Ex2 n -= 1 Ex -= x - K Ex2 -= ()x - K)#2def get_mean(): mundial K, n, Ex retorno K + Ex / ndef get_variance(): mundial n, Ex, Ex2 retorno ()Ex2 - Ex#2 / n) / ()n - 1)

Algoritmo de dos pasos

Un enfoque alternativo, que usa una fórmula diferente para la varianza, primero calcula la media de la muestra,

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

y luego calcula la suma de los cuadrados de las diferencias de la media,

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

donde s es la desviación estándar. Esto viene dado por el siguiente código:

def dos_pass_variance()datos): n = Len()datos) # = suma()datos) / n diferencia = suma[ ()x-#)#2 para x dentro datos ]) / ()n-1) retorno diferencia

Este algoritmo es numéricamente estable si n es pequeño. Sin embargo, los resultados de estos dos algoritmos simples ("ingenuo" y "dos pasos") pueden depender excesivamente del orden de los datos y pueden dar resultados deficientes para conjuntos de datos muy grandes. debido a errores repetidos de redondeo 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 un solo paso, inspeccionando cada valor $x_{i}$ sólo una vez; por ejemplo, cuando se recopilan los datos sin suficiente almacenamiento para mantener todos los valores, o cuando los costos del 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 las estadísticas requeridas se pueden calcular de forma numéricamente estable.

Las siguientes fórmulas se pueden utilizar para actualizar la diferencia media y (estimada) de la secuencia, para un elemento adicional x_n. Aquí, ${textstyle {overline {x}}_{n}={frac {1}{n}}sum _{i=1}^{n}x_{i}}$ denota la media de la muestra de la primera n muestras $(x_1,dots,x_n)$ , ${textstyle sigma _{n}^{2}={frac {1}{n}}sum _{i=1}^{n}left(x_{i}-{overline {x}}_{n}right)^{2}}$ su varianza de muestra sesgada, y ${textstyle s_{n}^{2}={frac {1}{n-1}}sum _{i=1}^{n}left(x_{i}-{overline {x}}_{n}right)^{2}}$ su varianza de muestra imparcial.

{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}}.}

1}" xmlns="http://www.w3.org/1998/Math/MathML">sn2=n− − 2n− − 1sn− − 12+()xn− − x̄ ̄ n− − 1)2n=sn− − 12+()xn− − x̄ ̄ n− − 1)2n− − sn− − 12n− − 1,n■1{displaystyle S_{n} {2}={frac {n-2}{n-1},s_{n-1}{2}+{frac {(x_{n}-{bar} {fn} {fn} {fn} {fn}fn} {fn}\fn}} {fn} {fn} {fn} {fn} {fn}fn} {fn}}}} {fn}}}} {fn} {fn} {fn}}}}} {fn}}}}}}}}}}}}}}}}}\\\n}}}}}}}}\\\\\\n}\\\\\\\\\\\\\\\\\\\fn1}}}}}\\\\\\\n}}}}}}}\\\\\\\\\\\\fn}}}}}}} {x}_{n-1} {n}} {n} {frac} {fn} {fn}} {fn}} {fn} {fn} {fn}} {fn} {fn} {fn}}} {fn} {fn}}}} {fn} {fn}}}}}}}} {f}}}}} {f}f}}}}}}}}}}}}}} {f} {f} {f}}}} {f}}}}}}}}}}}}}f} {f} {f}}f}}}}}}f} {fn} {f}fn}}}f}}}}}fn}f}f}}}}f}}}}}}}}}}}fn}}f}}}}}}}}f}}}}}}}}} {fn1}} {n-1}quad n]1}" aria-hidden="true" class="mwe-math-fallback-image-inline" src="https://wikimedia.org/api/rest_v1/media/math/render/svg/d3f407c573792f57f6a8b9799e70d9b13193bdc2" style="vertical-align: -2.005ex; width:75.186ex; height:6.176ex;"/>

Estas fórmulas sufren de inestabilidad numérica, ya que repetidamente restan un pequeño número de un gran número que escala con n. Una mejor cantidad para actualizar es la suma de cuadrados de diferencias de la media actual, ${textstyle sum _{i=1}^{n}(x_{i}-{bar {x}}_{n})^{2}}$ , aquí denotado $M_{2,n}$ :

{displaystyle {begin{aligned}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{aligned}}}

Este algoritmo fue encontrado por Welford, y se ha analizado a fondo. También es común denotar $M_k = bar x_k$ y $S_k = M_{2,k}$ .

A continuación se proporciona un ejemplo de implementación de Python para el algoritmo de Welford.

# Para un nuevo valor Valor, computar el nuevo recuento, nuevo medio, el nuevo M2.# significa acumular la media de todo el conjunto de datos# M2 agrega la distancia cuadrada de la media# Contar agrega el número de muestras vistas hasta ahoradef actualización()existentes Aggregate, nuevo Valor): ()Cuenta, #, M2) = existentes Aggregate Cuenta += 1 delta = nuevo Valor - # # += delta / Cuenta delta2 = nuevo Valor - # M2 += delta * delta2 retorno ()Cuenta, #, M2)# Recuperar la diferencia de media, diferencia y muestra de un agregadodef finalización()existentes Aggregate): ()Cuenta, #, M2) = existentes Aggregate si Cuenta . 2: retorno flotador()"nan") más: ()#, diferencia, muestra Diferencia) = ()#, M2 / Cuenta, M2 / ()Cuenta - 1) retorno ()#, diferencia, muestra Diferencia)

Este algoritmo es mucho menos propenso a la pérdida de 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 ciclo. 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 siguiente algoritmo paralelo ilustra cómo fusionar 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 los pesos vistos hasta ahora. West (1979) sugiere este algoritmo incremental:

def ponderado_incremental_varianza()data_weight_pairs): w_sum = w_sum2 = # = S = 0 para x, w dentro data_weight_pairs: w_sum = w_sum + w w_sum2 = w_sum2 + w#2 Quiero decir, viejo. = # # = Quiero decir, viejo. + ()w / w_sum) * ()x - Quiero decir, viejo.) S = S + w * ()x - Quiero decir, viejo.) * ()x - #) population_variance = S / w_sum # Corrección de Bessel para muestras ponderadas # Pesos de frecuencia sample_frequency_variance = S / ()w_sum - 1) # Pesos de fiabilidad sample_reliability_variance = S / ()w_sum - w_sum2 / w_sum)

Algoritmo paralelo

Chan et al. nota que el algoritmo en línea de Welford detallado arriba es un caso especial de un algoritmo que funciona para combinar conjuntos arbitrarios $A$ y $B$ :

{displaystyle {begin{aligned}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{aligned}}}

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 $n_A approx n_B$ y ambos son grandes, porque el error numérico en ${displaystyle delta ={bar {x}}_{B}-{bar {x}}_{A}}$ no se escala en la forma en que está en $n_B = 1$ caso. En tales casos, prefiera ${textstyle {bar {x}}_{AB}={frac {n_{A}{bar {x}}_{A}+n_{B}{bar {x}}_{B}}{n_{AB}}}}$ .

def paralel_variancia()n, avg_a, M2_a, n_b, avg_b, M2_b): n = n + n_b delta = avg_b - avg_a M2 = M2_a + M2_b + delta # 2 * n * n_b / n var_ab = M2 / ()n - 1) retorno var_ab

Esto se puede generalizar para permitir la paralelización con AVX, con GPU y clústeres de computadoras, y para la covarianza.

Ejemplo

Suponga que todas las operaciones de punto flotante utilizan aritmética de doble precisión estándar IEEE 754. Considere la muestra (4, 7, 13, 16) de una población infinita. Según esta muestra, la media poblacional estimada es 10 y la estimación imparcial 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⁸ + 4, 10⁸ + 7, 10⁸ + 13, 10⁸ + 16), lo 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,333333333333332 en lugar de 30.

Si bien esta pérdida de precisión puede ser tolerable y vista como una falla menor del algoritmo ingenuo, aumentar aún más la compensación hace que el error sea catastrófico. Considere la muestra (10⁹ + 4, 10⁹ + 7, 10⁹ + 13, 10⁹ + 16). Una vez más, la varianza estimada de la población 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 serio 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 amplía las fórmulas de Chan para calcular los momentos centrales tercero y cuarto, necesarios, por ejemplo, al estimar la asimetría y la curtosis:

{displaystyle {begin{aligned}M_{3,X}=M_{3,A}+M_{3,B}&{}+delta ^{3}{frac {n_{A}n_{B}(n_{A}-n_{B})}{n_{X}^{2}}}+3delta {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}left(n_{A}^{2}-n_{A}n_{B}+n_{B}^{2}right)}{n_{X}^{3}}}\[6pt]&{}+6delta ^{2}{frac {n_{A}^{2}M_{2,B}+n_{B}^{2}M_{2,A}}{n_{X}^{2}}}+4delta {frac {n_{A}M_{3,B}-n_{B}M_{3,A}}{n_{X}}}end{aligned}}}

Aquí está. $M_k$ son otra vez las sumas de poderes de diferencias de la media ${textstyle sum (x-{overline {x}})^{k}}$ , dar

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

Para el caso incremental (es decir, $B = {x}$ ), esto simplifica:

{displaystyle {begin{aligned}delta &=x-m\[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 {3delta M_{2}}{n}}\[5pt]M_{4}'&=M_{4}+{frac {delta ^{4}(n-1)(n^{2}-3n+3)}{n^{3}}}+{frac {6delta ^{2}M_{2}}{n^{2}}}-{frac {4delta M_{3}}{n}}end{aligned}}}

Al preservar el valor $delta / n$ , sólo se necesita una operación de división y las estadísticas de orden superior se pueden calcular por lo tanto para un pequeño costo incremental.

Un ejemplo del algoritmo en línea para la curtosis implementado como se describe es:

def online_kurtosis()datos): n = # = M2 = M3 = M4 = 0 para x dentro datos: n1 = n n = n + 1 delta = x - # delta_n = delta / n delta_n2 = delta_n # 2 mandato1 = delta * delta_n * n1 # = # + delta_n M4 = M4 + mandato1 * delta_n2 * ()n#2 - 3*n + 3) + 6 * delta_n2 * M2 - 4 * delta_n * M3 M3 = M3 + mandato1 * delta_n * ()n - 2) - 3 * delta_n * M2 M2 = M2 + mandato1 # Nota, también puede calcular la varianza utilizando M2, y la asiduidad usando M3 # Precaución: Si todas las entradas son iguales, M2 será 0, resultando en una división por 0. kurtosis = ()n * M4) / ()M2 # 2) - 3 retorno kurtosis

Pébaÿ extiende aún más estos resultados a momentos centrales de orden arbitrario, para los casos incremental y por pares, y posteriormente Pébaÿ et al. para momentos ponderados y compuestos. También se pueden encontrar allí fórmulas similares para la covarianza.

Choi y Sweetman ofrecen dos métodos alternativos para calcular la asimetría y la curtosis, cada uno de los cuales puede ahorrar requisitos sustanciales de memoria de 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 se convierte efectivamente en un algoritmo de un solo paso para momentos más altos. Un beneficio es que los cálculos de momentos estadísticos pueden llevarse a cabo con una precisión arbitraria, de modo que los cálculos pueden ajustarse a 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 representa gráficamente 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:

H(x_k)=frac{h(x_k)}{A}

Donde $h(x_k)$ y $H(x_k)$ representan la frecuencia y la frecuencia relativa en bin $x_{k}$ y ${textstyle A=sum _{k=1}^{K}h(x_{k}),Delta x_{k}}$ es el área total del histograma. Después de esta normalización, la $n$ momentos crudos y momentos centrales $x(t)$ se puede calcular a partir del histograma relativo:

{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)}{Big)}^{n},H(x_{k}),Delta x_{k}={frac {1}{A}}sum _{k=1}^{K}{Big (}x_{k}-m_{1}^{(h)}{Big)}^{n}h(x_{k}),Delta x_{k}}

donde el superscript $^{(h)}$ indica que los momentos se calculan a partir del histograma. Para el ancho constante del bin ${displaystyle Delta x_{k}=Delta x}$ estas dos expresiones se pueden simplificar usando $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)}{Big)}^{n}h(x_{k})}

El segundo enfoque de Choi y Sweetman es una metodología analítica para combinar momentos estadísticos de segmentos individuales de una historia temporal de modo que los momentos generales resultantes sean los de la historia temporal completa. Esta metodología podría usarse 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 $Q$ conjuntos de momentos estadísticos son conocidos: $(gamma_{0,q},mu_{q},sigma^2_{q},alpha_{3,q},alpha_{4,q}) quad$ para ${displaystyle q=1,2,ldotsQ}$ , entonces cada uno $gamma _{n}$ puede se expresa en términos del equivalente $n$ momentos crudos:

{displaystyle gamma _{n,q}=m_{n,q}gamma _{0,q}qquad quad {textrm {for}}quad n=1,2,3,4quad {text{ and }}quad q=1,2,dotsQ}

Donde $gamma_{0,q}$ se toma generalmente como la duración de la $q^{th}$ tiempo-historia, o el número de puntos si $Delta t$ es constante.

El beneficio de expresar los momentos estadísticos en términos de $gamma$ es que $Q$ conjuntos se pueden combinar por adición, y no hay límite superior en el valor de $Q$ .

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

donde el subscript $_c$ representa la historia del tiempo concatenado o combinado $gamma$ . Estos valores combinados $gamma$ puede entonces ser transformado inversamente en momentos crudos que representan la historia del tiempo completo

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

Relaciones conocidas entre los momentos crudos $m_{n}$ ) y los momentos centrales ( ${displaystyle theta _{n}=operatorname {E} [(x-mu)^{n}])}$ ) son utilizados para calcular los momentos centrales de la historia del tiempo concatenado. Finalmente, los momentos estadísticos de la historia concatenada se computan desde los momentos centrales:

{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}

Covarianza

Se pueden usar algoritmos muy similares para calcular la covarianza.

Algoritmo ingenuo

El algoritmo ingenuo es

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

Para el algoritmo anterior, se podría usar el siguiente código de Python:

def inive_covariance()data1, data2): n = Len()data1) suma1 = suma()data1) sum2 = suma()data2) sum12 = suma[ i1*i2 para i1,i2 dentro Cierre()data1,data2) ]) covariancia = ()sum12 - suma1 * sum2 / n) / n retorno covariancia

Con estimación de la media

En cuanto a la varianza, la covariancia de dos variables aleatorias también es invariable de cambio, así que dadas dos valores constantes $k_x$ y ${displaystyle k_{y},}$ puede ser escrito:

{displaystyle operatorname {Cov} (X,Y)=operatorname {Cov} (X-k_{x},Y-k_{y})={dfrac {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}}.}

y nuevamente elegir un valor dentro del rango de valores estabilizará la fórmula contra cancelaciones catastróficas y la hará más sólida contra grandes sumas. Tomando el primer valor de cada conjunto de datos, el algoritmo se puede escribir como:

def cambio_data_covariancia()data_x, data_y): n = Len()data_x) si n . 2: retorno 0 kx = data_x[0] ky = data_y[0] Ex = Ey = Exy = 0 para ix, i dentro Cierre()data_x, data_y): Ex += ix - kx Ey += i - ky Exy += ()ix - kx) * ()i - ky) retorno ()Exy - Ex * Ey / n) / n

Dos pases

El algoritmo de dos pasos primero calcula las medias de la muestra y luego la covarianza:

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

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

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

El algoritmo de dos pasos se puede escribir como:

def dos_pass_covariance()data1, data2): n = Len()data1) media1 = suma()data1) / n media2 = suma()data2) / n covariancia = 0 para i1, i2 dentro Cierre()data1, data2): a = i1 - media1 b = i2 - media2 covariancia += a * b / n retorno covariancia

Una versión compensada ligeramente más precisa realiza el algoritmo ingenuo completo en los residuos. Las sumas finales ${textstyle sum_i x_i}$ y ${textstyle sum _{i}y_{i}}$ debería ser cero, pero el segundo paso compensa cualquier error pequeño.

En línea

Existe un algoritmo estable de un paso, similar al algoritmo en línea para calcular la varianza, que calcula el co-momento ${textstyle C_{n}=sum _{i=1}^{n}(x_{i}-{bar {x}}_{n})(y_{i}-{bar {y}}_{n})}$ :

{displaystyle {begin{alignedat}{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{alignedat}}}

La asimetría aparente en esa última ecuación se debe al hecho de que ${textstyle (x_{n}-{bar {x}}_{n})={frac {n-1}{n}}(x_{n}-{bar {x}}_{n-1})}$ , por lo que ambos términos de actualización son iguales ${textstyle {frac {n-1}{n}}(x_{n}-{bar {x}}_{n-1})(y_{n}-{bar {y}}_{n-1})}$ . Incluso mayor precisión se puede lograr por primera vez computar los medios, luego utilizando el algoritmo estable de un paso en los residuos.

Por lo tanto, la covarianza se puede calcular como

{displaystyle {begin{aligned}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{aligned}}}

def online_covariancia()data1, data2): # = mezquino = C = n = 0 para x, Sí. dentro Cierre()data1, data2): n += 1 dx = x - # # += dx / n mezquino += ()Sí. - mezquino) / n C += dx * ()Sí. - mezquino) population_covar = C / n # Corrección de Bessel para la varianza de muestra sample_covar = C / ()n - 1)

También se puede hacer una pequeña modificación para calcular la covarianza ponderada:

def online_pesado_covariancia()data1, data2, data3): # = mezquino = 0 wsum = wsum2 = 0 C = 0 para x, Sí., w dentro Cierre()data1, data2, data3): wsum += w wsum2 += w * w dx = x - # # += ()w / wsum) * dx mezquino += ()w / wsum) * ()Sí. - mezquino) C += w * dx * ()Sí. - mezquino) population_covar = C / wsum # Corrección de Bessel para la varianza de muestra # Pesos de frecuencia sample_frequency_covar = C / ()wsum - 1) # Pesos de fiabilidad sample_reliability_covar = C / ()wsum - wsum2 / wsum)

Del mismo modo, existe una fórmula para combinar las covarianzas de dos conjuntos que se puede usar para paralelizar el 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}}}.}

Versión por lotes ponderada

Una versión del algoritmo en línea ponderado que hace batched actualizado también existe: dejar ${displaystyle w_{1},dots w_{N}}$ denota los pesos, y escribe