Algoritmos para calcular la varianza

Los algoritmos para calcular la varianza desempeñan 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 una inestabilidad numérica , así como a un desbordamiento aritmético cuando se trabaja con valores grandes.

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

Utilizando 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 el siguiente:

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 hacer 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 lo 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.

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

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

con cualquier constante, lo que conduce a la nueva fórmula ${\estilo de visualización K}$

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

Cuanto más cerca esté del valor medio, más preciso será el resultado, pero simplemente eligiendo un valor dentro del rango de las muestras se garantizará la estabilidad deseada. Si los valores son pequeños, no hay problemas con la suma de sus cuadrados; por el contrario, si son grandes, significa necesariamente 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]${\estilo de visualización K}$${\displaystyle (x_{i}-K)}$

Si solo se toma la primera muestra, el algoritmo se puede escribir en lenguaje de programación Python como ${\estilo de visualización 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  +=  1 
        Ex  +=  x  -  K 
        Ex2  +=  ( x  -  K )  **  2 
    variance  =  ( 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) si los datos son muestras de una población más grande 
    devuelve  varianza

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

K  =  Ex  =  Ex2  =  0,0 
n  =  0


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

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

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

def  obtener_varianza (): 
    global  n ,  Ex ,  Ex2 
    devuelve  ( Ex2  -  Ex ** 2  /  n )  /  ( n  -  1 )

Un enfoque alternativo, que utiliza 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 con respecto a la media,

${\displaystyle {\text{varianza de muestra}}=s^{2}={\dfrac {\sum _{i=1}^{n}(x_{i}-{\bar {x}})^{2}}{n-1}},}$

donde s es la desviación estándar. Esto se obtiene mediante el siguiente código:

def  two_pass_variance ( datos ): 
    n  =  len ( datos ) 
    media  =  suma ( datos )  /  n 
    varianza  =  suma (( x  -  media )  **  2  para  x  en  datos )  /  ( n  -  1 ) 
    devuelve  varianza

Este algoritmo es numéricamente estable si n es pequeño. ^{[1] [4]} Sin embargo, los resultados de ambos 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 error de redondeo repetido en la acumulación de las sumas. Se pueden utilizar técnicas como la suma compensada para combatir este error hasta cierto punto.

A menudo resulta útil poder calcular la varianza en una sola pasada , inspeccionando cada valor solo una vez; por ejemplo, cuando se recopilan los datos sin suficiente almacenamiento para guardar todos los valores, o cuando los costos de acceso a la memoria dominan los de cálculo. Para un algoritmo en línea de este tipo , se requiere una relación de recurrencia entre las cantidades a partir de las cuales se pueden calcular las estadísticas requeridas de una manera numéricamente estable. $Estilo de visualización x_{i}}$

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í, denota la media muestral de las primeras n muestras , su varianza muestral sesgada y su varianza muestral no sesgada . ${\textstyle {\overline {x}}_{n}={\frac {1}{n}}\sum _{i=1}^{n}x_{i}}$${\displaystyle (x_{1},\puntos ,x_{n})}$${\textstyle \sigma _{n}^{2}={\frac {1}{n}}\sum _{i=1}^{n}\left(x_{i}-{\overline {x}}_{n}\right)^{2}}$${\textstyle 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}$

Estas fórmulas sufren inestabilidad numérica ^{[ cita requerida ]} , ya que restan repetidamente un número pequeño de un número grande que escala con n . Una cantidad mejor para actualizar es la suma de los cuadrados de las diferencias con respecto a la media actual, , que aquí se denota como : ${\textstyle \suma _{i=1}^{n}(x_{i}-{\bar {x}}_{n})^{2}}$$Estilo de visualización 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 descubierto por Welford, ^{[5] [6]} y ha sido analizado exhaustivamente. ^{[2] [7]} También es común denotar y . ^[8]${\displaystyle M_{k}={\bar {x}}_{k}}$${\displaystyle S_{k}=M_{2,k}}$

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

# Para un nuevo valor new_value, calcula el nuevo conteo, la nueva media, el nuevo M2. 
# media acumula la media de todo el conjunto de datos 
# M2 agrega la distancia al cuadrado desde la media 
# count agrega la cantidad de muestras vistas hasta el momento 
def  update ( existing_gregate ,  new_value ): 
    ( count ,  mean ,  M2 )  =  existing_gregate 
    count  +=  1 
    delta  =  new_value  -  mean 
    mean  +=  delta  /  count 
    delta2  =  new_value  -  mean 
    M2  +=  delta  *  delta2 
    return  ( count ,  mean ,  M2 )

# Recuperar la media, la varianza y la varianza de la muestra de un agregado 
def  finalize ( existing_gregate ): 
    ( count ,  mean ,  M2 )  =  existing_gregate 
    if  count  <  2 : 
        return  float ( "nan" ) 
    else : 
        ( mean ,  variance ,  sample_variance )  =  ( mean ,  M2  /  count ,  M2  /  ( count  -  1 )) 
        return  ( mean ,  variance ,  sample_variance )

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 bucle. Para un algoritmo de dos pasadas 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 fusionar múltiples conjuntos de estadísticas calculadas en línea.

El algoritmo se puede ampliar para manejar pesos de muestra desiguales, reemplazando el simple contador n con la suma de pesos vistos hasta ahora. West (1979) ^[9] sugiere este algoritmo incremental :

def  varianza_incremental_ponderada ( pares_ponderados_de_datos ): 
    suma_w  =  suma_w2  =  media  =  S  =  0

    para  x ,  w  en  pares_de_pesos_de_datos : 
        w_suma  =  w_suma  +  w 
        w_suma2  =  w_suma2  +  w ** 2 
        media_antigua  =  media media 
        = media_antigua + ( w / w_suma ) * ( x - media_antigua ) S = S + w * ( x - media_antigua ) * ( x - media )          
                    

    population_variance  =  S  /  w_sum 
    # Corrección de Bessel para muestras ponderadas 
    # para ponderaciones de frecuencia entera 
    sample_frequency_variance  =  S  /  ( w_sum  -  1 )

Chan et al. ^[10] señalan que el algoritmo en línea de Welford detallado anteriormente es un caso especial de un algoritmo que funciona para combinar conjuntos arbitrarios y : ${\estilo de visualización A}$${\estilo de visualización 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 y ambos son grandes, porque el error numérico en no se reduce de la forma en que lo hace en este caso. En tales casos, es preferible . ${\displaystyle n_{A}\approx n_{B}}$${\displaystyle \delta ={\bar {x}}_{B}-{\bar {x}}_{A}}$${\displaystyle n_{B}=1}$${\textstyle {\bar {x}}_{AB}={\frac {n_{A}{\bar {x}}_{A}+n_{B}{\bar {x}}_{B}}{n_{AB}}}}$

def  varianza_paralela ( 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 ) 
    devolver  var_ab

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

Supongamos que todas las operaciones de punto flotante utilizan la aritmética de doble precisión estándar IEEE 754. Consideremos 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 no sesgada de la varianza poblacional es 30. Tanto el algoritmo ingenuo como el algoritmo de dos pasadas calculan estos valores correctamente.

A continuación, considere la muestra ( 10 ⁸  + 4 , 10 ⁸  + 7 , 10 ⁸  + 13 , 10 ⁸  + 16 ), que da lugar a la misma varianza estimada que la primera muestra. El algoritmo de dos pasadas 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 verse como un defecto menor del algoritmo ingenuo, aumentar aún más el desfase hace que el error sea catastrófico. Considere la muestra ( 10 ⁹  + 4 , 10 ⁹  + 7 , 10 ⁹  + 13 , 10 ⁹  + 16 ). Nuevamente, la varianza poblacional estimada de 30 se calcula correctamente mediante el algoritmo de dos pasadas, pero el algoritmo ingenuo ahora la calcula como −170,66666666666666. Este es un problema grave 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.

Terriberry ^[11] extiende las fórmulas de Chan para calcular el tercer y cuarto momento central , 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}}}+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}\left(n_{A}^{2}-n_{A}n_{B}+n_{B}^{2}\right)}{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{aligned}}}$

Aquí están nuevamente las sumas de potencias de diferencias de la media , dando ${\displaystyle M_{k}}$${\textstyle \sum (x-{\overline {x}})^{k}}$

${\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, ), esto se simplifica a: ${\displaystyle B=\{x\}}$

${\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 {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{aligned}}}$

Al preservar el valor , solo se necesita una operación de división y las estadísticas de orden superior se pueden calcular con un costo incremental mínimo. ${\displaystyle \delta /n}$

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

def  online_kurtosis ( datos ): 
    n  =  media  =  M2  =  M3  =  M4  =  0

    para  x  en  los datos : 
        n1  =  n 
        n  =  n  +  1 
        delta  =  x  -  media 
        delta_n  =  delta  /  n 
        delta_n2  =  delta_n ** 2 
        término1  =  delta  *  delta_n  *  n1 
        media  =  media  +  delta_n 
        M4  =  M4  +  término1  *  delta_n2  *  ( n ** 2  -  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

    # Nota, 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 por 0. 
    kurtosis  =  ( n  *  M4 )  /  ( M2 ** 2 )  -  3 
    return  kurtosis

Pébaÿ ^{[12] 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. ^[13] para momentos ponderados y compuestos. También se pueden encontrar allí 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 requisitos sustanciales de memoria de computadora y tiempo de CPU en ciertas aplicaciones. El primer enfoque es calcular los momentos estadísticos separando los datos en contenedores y luego calculando los momentos a partir de la geometría del histograma resultante, lo que efectivamente se convierte en un algoritmo de una sola pasada para momentos más altos. Un beneficio es que los cálculos de momentos estadísticos se pueden realizar con una precisión arbitraria de modo que los cálculos se puedan ajustar a la precisión de, por ejemplo, el formato de almacenamiento de datos o el hardware de medición original. Un histograma relativo de una variable aleatoria se puede construir de la manera convencional: el rango de valores potenciales se divide en contenedores y se cuenta y se grafica el número de ocurrencias dentro de cada contenedor de modo que el área de cada rectángulo sea igual a la porción de los valores de muestra dentro de ese contenedor:

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

donde y representan la frecuencia y la frecuencia relativa en el intervalo y es el área total del histograma. Después de esta normalización, los momentos brutos y los momentos centrales de se pueden calcular a partir del histograma relativo: ${\displaystyle h(x_{k})}$${\displaystyle H(x_{k})}$${\displaystyle x_{k}}$${\textstyle 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)}{\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 superíndice indica que los momentos se calculan a partir del histograma. Para un ancho de intervalo constante, estas dos expresiones se pueden simplificar utilizando : ${\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)}{\Big )}^{n}h(x_{k})}$

El segundo enfoque de Choi y Sweetman ^[14] es una metodología analítica para combinar momentos estadísticos de segmentos individuales de un historial temporal de modo que los momentos generales resultantes sean los del historial temporal completo. Esta metodología podría utilizarse para el cálculo paralelo de momentos estadísticos con la posterior combinación de esos momentos, o para la combinación de momentos estadísticos calculados en momentos secuenciales.

Si se conocen conjuntos de momentos estadísticos: para , entonces cada uno puede expresarse en términos de los 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 {for}}\quad n=1,2,3,4\quad {\text{ and }}\quad q=1,2,\dots ,Q}$

donde generalmente se toma como la duración del historial de tiempo o el número de puntos si es constante. ${\displaystyle \gamma _{0,q}}$${\displaystyle q^{th}}$${\displaystyle \Delta t}$

La ventaja de expresar los momentos estadísticos en términos de es que los conjuntos se pueden combinar mediante adición y no hay un límite superior en el 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}$

donde el subíndice representa el historial temporal concatenado o combinado . Estos valores combinados de pueden luego transformarse inversamente en momentos sin procesar que representan el historial temporal concatenado completo ${\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}$

Las relaciones conocidas entre los momentos brutos ( ) y los momentos centrales ( ) se utilizan luego para calcular los momentos centrales de la historia temporal concatenada. Finalmente, los momentos estadísticos de la historia concatenada se calculan a partir de los momentos centrales: ${\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}$

Se pueden utilizar algoritmos muy similares para calcular la covarianza .

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 utilizar el siguiente código Python:

def  naive_covariance ( datos1 ,  datos2 ): 
    n  =  len ( datos1 ) 
    suma1  =  suma ( datos1 ) 
    suma2  =  suma ( datos2 ) 
    suma12  =  suma ([ i1  *  i2  para  i1 ,  i2  en  zip ( datos1 ,  datos2 )])

    covarianza  =  ( suma12  -  suma1  *  suma2  /  n )  /  n 
    devuelve  covarianza

En cuanto a la varianza, la covarianza de dos variables aleatorias también es invariante al desplazamiento, por lo que dados dos valores constantes cualesquiera , se puede escribir: ${\displaystyle k_{x}}$${\displaystyle k_{y},}$

${\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, de nuevo, si se elige un valor dentro del rango de valores, se estabilizará la fórmula frente a una cancelación catastrófica y se la hará más robusta frente a grandes sumas. Si se toma el primer valor de cada conjunto de datos, el algoritmo se puede escribir de la siguiente manera:

def  shifted_data_covariance ( datos_x ,  datos_y ): 
    n  =  len ( datos_x ) 
    si  n  <  2 : 
        devuelve  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 ) 
    devuelve  ( Exy  -  Ex  *  Ey  /  n )  /  n

El algoritmo de dos pasos primero calcula las medias de 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  covarianza_de_dos_pasos ( datos1 ,  datos2 ): 
    n  =  len ( datos1 ) 
    media1  =  suma ( datos1 )  /  n 
    media2  =  suma ( datos2 )  /  n

    covarianza  =  0 
    para  i1 ,  i2  en  zip ( datos1 ,  datos2 ): 
        a  =  i1  -  media1 
        b  =  i2  -  media2 
        covarianza  +=  a  *  b  /  n 
    devuelve  covarianza

Una versión compensada ligeramente más precisa ejecuta el algoritmo ingenuo completo sobre los residuos. Las sumas finales deberían ser cero, pero la segunda pasada compensa cualquier pequeño error. ${\textstyle \sum _{i}x_{i}}$${\textstyle \sum _{i}y_{i}}$

Existe un algoritmo estable de una sola pasada, 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 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 las medias y luego utilizando el algoritmo estable de una sola pasada en los residuos. ${\textstyle (x_{n}-{\bar {x}}_{n})={\frac {n-1}{n}}(x_{n}-{\bar {x}}_{n-1})}$${\textstyle {\frac {n-1}{n}}(x_{n}-{\bar {x}}_{n-1})(y_{n}-{\bar {y}}_{n-1})}$

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_covariance ( datos1 ,  datos2 ): 
    mediax  =  mediay  =  C  =  n  =  0 
    para  x ,  y  en  zip ( datos1 ,  datos2 ): 
        n  +=  1 
        dx  =  x  -  mediax 
        mediax  +=  dx  /  n 
        mediay  +=  ( y  -  mediay )  /  n 
        C  +=  dx  *  ( y  -  mediay )

    población_covar  =  C  /  n 
    # Corrección de Bessel para la varianza de la muestra 
    muestra_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 ,  datos3 ): 
    mediax  =  mediay  =  0 
    wsum  =  wsum2  =  0 
    C  =  0 
    para  x ,  y ,  w  en  zip ( datos1 ,  datos2 ,  datos3 ): 
        wsum  +=  w 
        wsum2  +=  w  *  w 
        dx  =  x  -  mediax 
        mediax  +=  ( w  /  wsum )  *  dx 
        mediay  +=  ( w  /  wsum )  *  ( y  -  mediay ) 
        C  +=  w  *  dx  *  ( y  -  mediay )

    population_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 puede utilizarse para paralelizar el cálculo: ^[3]

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

También existe una versión del algoritmo ponderado en línea que realiza actualizaciones por lotes: denotemos los pesos y escribamos ${\displaystyle w_{1},\dots w_{N}}$

${\displaystyle {\begin{alignedat}{2}{\bar {x}}_{n+k}&={\bar {x}}_{n}&\,+\,&{\frac {\sum _{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{alignedat}}}$

La covarianza puede entonces calcularse como

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

• Algoritmo de suma de Kahan
• Desviaciones al cuadrado de la media
• Método Yamartino

• Weisstein, Eric W. "Cálculo de la varianza de la muestra". MathWorld .