Algoritmos para calcular la varianza - Algorithms for calculating variance

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:

${\ Displaystyle \ sigma ^ {2} = {\ overline {(x ^ {2})}} - {\ bar {x}} ^ {2} = \ displaystyle {\ 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 insesgada 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 {\ suma _ {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:

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

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.

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

con cualquier constante, lo que conduce a la nueva fórmula ${\ Displaystyle K}$

${\ Displaystyle \ sigma ^ {2} = \ displaystyle {\ 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 elegir un valor dentro del rango de muestras 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 de la fórmula es siempre menor que el primero, por lo que no puede producirse ninguna cancelación. ${\ Displaystyle K}$${\ Displaystyle (x_ {i} -K)}$

Si solo se toma la primera muestra, ya que el algoritmo se puede escribir en el lenguaje de programación Python como ${\ Displaystyle 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 = 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 también facilita el cálculo incremental que se puede expresar como

K = n = Ex = Ex2 = 0.0

def add_variable(x):
    global K, n, Ex, Ex2
    if 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 / n

def get_variance():
    global n, 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,

${\ 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 {varianza de muestra}} = s ^ {2} = \ displaystyle {\ frac {\ 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 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 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 datos. 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 valor solo 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. ${\ Displaystyle 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 insesgada . ${\ Displaystyle {\ overline {x}} _ {n} = {\ frac {1} {n}} \ sum _ {i = 1} ^ {n} x_ {i}}$${\ Displaystyle (x_ {1}, \ dots, x_ {n})}$${\ Displaystyle \ sigma _ {n} ^ {2} = {\ frac {1} {n}} \ sum _ {i = 1} ^ {n} \ left (x_ {i} - {\ overline {x} } _ {n} \ derecha) ^ {2}}$${\ Displaystyle s_ {n} ^ {2} = {\ frac {1} {n-1}} \ sum _ {i = 1} ^ {n} \ left (x_ {i} - {\ overline {x} } _ {n} \ derecha) ^ {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} }{norte}}.}$

${\ 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 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í denotada : ${\ Displaystyle \ textstyle \ sum _ {i = 1} ^ {n} (x_ {i} - {\ bar {x}} _ {n}) ^ {2}}$${\ Displaystyle M_ {2, n}}$

${\ Displaystyle {\ begin {alineado} 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 {alineado}}}$

Este algoritmo fue encontrado por Welford y ha sido analizado a fondo. También es común denotar y . ${\ Displaystyle M_ {k} = {\ bar {x}} _ {k}}$${\ Displaystyle S_ {k} = M_ {2, k}}$

A continuación se proporciona un ejemplo de implementación de 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
    if count < 2:
        return float("nan")
    else:
        (mean, variance, sampleVariance) = (mean, M2 / count, M2 / (count - 1))
        return (mean, variance, sampleVariance)

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) sugiere este algoritmo incremental :

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 * w
        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
    # Bessel's correction for weighted samples
    # Frequency weights
    sample_frequency_variance = S / (w_sum - 1)
    # Reliability weights
    sample_reliability_variance = S / (w_sum - w_sum2 / w_sum)

Algoritmo paralelo

Chan y col. 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 : ${\ Displaystyle A}$${\ Displaystyle B}$

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

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 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 . ${\ Displaystyle n_ {A} \ approx n_ {B}}$${\ Displaystyle \ delta = {\ bar {x}} _ {B} - {\ bar {x}} _ {A}}$${\ Displaystyle n_ {B} = 1}$${\ Displaystyle {\ bar {x}} _ {AB} = {\ frac {n_ {A} {\ bar {x}} _ {A} + n_ {B} {\ bar {x}} _ {B} }{coger}}}}$

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.

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 ⁸  + 4 , 10 ⁸  + 7 , 10 ⁸  + 13 , 10 ⁸  + 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 ⁹  + 4 , 10 ⁹  + 7 , 10 ⁹  + 13 , 10 ⁹  + 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 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 :

${\ Displaystyle {\ begin {alineado} 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} \ izquierda (n_ {A} ^ {2} -n_ {A} n_ {B} + n_ {B} ^ {2} \ derecha)} {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 {alineado}}}$

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

${\ Displaystyle {\ begin {alineado} & {\ text {sesgo}} = g_ {1} = {\ frac {{\ sqrt {n}} M_ {3}} {M_ {2} ^ {3/2} }}, \\ [4pt] & {\ text {curtosis}} = g_ {2} = {\ frac {nM_ {4}} {M_ {2} ^ {2}}} - 3. \ end {alineado} }}$

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

${\ Displaystyle {\ begin {alineado} \ delta & = xm \\ [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 {alineado}}}$

Al conservar 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. ${\ Displaystyle \ delta / n}$

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

    # Note, you may also calculate variance using M2, and skewness using M3
    # Caution: If all the inputs are the same, M2 will be 0, resulting in a division by 0.
    kurtosis = (n * M4) / (M2 * M2) - 3
    return kurtosis

Pébaÿ extiende estos resultados a momentos centrales de orden arbitrario , para los casos incrementales y por pares, y posteriormente Pébaÿ et al. para momentos ponderados y compuestos. También se pueden encontrar 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 importantes requisitos de memoria de la 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, que efectivamente se convierte en un algoritmo de una pasada para momentos superiores. Un beneficio es que los cálculos de momento estadístico 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:

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

donde y representan la frecuencia y la frecuencia relativa en bin y es el área total del histograma. Después de esta normalización, los momentos sin procesar y los momentos centrales de se pueden calcular a partir del histograma relativo: ${\ Displaystyle h (x_ {k})}$${\ Displaystyle H (x_ {k})}$${\ Displaystyle x_ {k}}$${\ Displaystyle 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)} {\ Grande)} ^ {n} \, H (x_ {k}) \, \ Delta x_ {k} = {\ frac {1} {A}} \ sum _ {k = 1} ^ {K} {\ Grande (} 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 contenedor constante, estas dos expresiones se pueden simplificar usando : ${\ 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 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: 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 {para}} \ quad n = 1,2,3,4 \ quad {\ text {y}} \ 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}$

El beneficio de expresar los momentos estadísticos en términos de es que los conjuntos se pueden combinar mediante la suma y no existe 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 {para}} n = 0,1,2 , 3,4}$

donde el subíndice representa la historia de tiempo concatenada o combinada . Estos valores combinados de pueden entonces transformarse inversamente en momentos brutos que representan la historia del tiempo concatenada completa. ${\ Displaystyle _ {c}}$${\ Displaystyle \ gamma}$${\ Displaystyle \ gamma}$

${\ Displaystyle m_ {n, c} = {\ frac {\ gamma _ {n, c}} {\ gamma _ {0, c}}} \ quad {\ text {para}} 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 de tiempo 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}$

Covarianza

Se pueden utilizar algoritmos muy similares para calcular la covarianza .

Algoritmo ingenuo

El algoritmo ingenuo es:

${\ Displaystyle \ operatorname {Cov} (X, Y) = \ displaystyle {\ 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 naive_covariance(data1, data2):
    n = len(data1)
    sum12 = 0
    sum1 = sum(data1)
    sum2 = sum(data2)

    for i1, i2 in zip(data1, data2):
        sum12 += i1 * i2

    covariance = (sum12 - sum1 * sum2 / n) / n
    return covariance

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: ${\ Displaystyle k_ {x}}$${\ Displaystyle k_ {y},}$

${\ Displaystyle \ operatorname {Cov} (X, Y) = \ operatorname {Cov} (X-k_ {x}, Y-k_ {y}) = \ displaystyle {\ frac {\ 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, 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 shifted_data_covariance(data_x, data_y):
    n = len(data_x)
    if n < 2:
        return 0
    kx = data_x[0]
    ky = data_y[0]
    Ex = Ey = Exy = 0
    for ix, iy in zip(data_x, data_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:

${\ Displaystyle {\ bar {x}} = \ Displaystyle \ sum _ {i = 1} ^ {n} x_ {i} / n}$

${\ Displaystyle {\ bar {y}} = \ Displaystyle \ sum _ {i = 1} ^ {n} y_ {i} / n}$

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

El algoritmo de dos pasos se puede escribir como:

def two_pass_covariance(data1, data2):
    n = len(data1)

    mean1 = sum(data1) / n
    mean2 = sum(data2) / n

    covariance = 0

    for i1, i2 in 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 deben ser cero, pero la segunda pasada compensa cualquier pequeño error. ${\ Displaystyle \ textstyle \ sum x_ {i}}$${\ Displaystyle \ textstyle \ sum y_ {i}}$

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 : ${\ Displaystyle \ textstyle C_ {n} = \ sum _ {i = 1} ^ {n} (x_ {i} - {\ bar {x}} _ {n}) (y_ {i} - {\ bar { y}} _ {n})}$

${\ displaystyle {\ begin {alineado} {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 {alineado}}}$

La aparente asimetría en esa última ecuación se debe al hecho de que , por lo tanto, 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. ${\ Displaystyle \ textstyle (x_ {n} - {\ bar {x}} _ {n}) = {\ frac {n-1} {n}} (x_ {n} - {\ bar {x}} _ {n-1})}$${\ Displaystyle \ textstyle {\ frac {n-1} {n}} (x_ {n} - {\ bar {x}} _ {n-1}) (y_ {n} - {\ bar {y}} _ {n-1})}$

Por tanto, la covarianza se puede calcular como

${\ Displaystyle {\ begin {alineado} \ 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 {alineado}}}$

def online_covariance(data1, data2):
    meanx = meany = C = n = 0
    for x, y in zip(data1, data2):
        n += 1
        dx = x - meanx
        meanx += dx / n
        meany += (y - meany) / n
        C += dx * (y - meany)

    population_covar = C / n
    # Bessel's correction for sample variance
    sample_covar = C / (n - 1)

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

def online_weighted_covariance(data1, data2, data3):
    meanx = meany = 0
    wsum = wsum2 = 0
    C = 0
    for x, y, w in zip(data1, data2, data3):
        wsum += w
        wsum2 += w * w
        dx = x - meanx
        meanx += (w / wsum) * dx
        meany += (w / wsum) * (y - meany)
        C += w * dx * (y - meany)

    population_covar = C / wsum
    # Bessel's correction for sample variance
    # Frequency weights
    sample_frequency_covar = C / (wsum - 1)
    # Reliability weights
    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:

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

También existe una versión del algoritmo ponderado en línea que se actualiza por lotes: denote los pesos y escriba ${\ Displaystyle w_ {1}, \ dots w_ {N}}$

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

La covarianza se puede calcular como

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

Ver también

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

Referencias

enlaces externos

• Weisstein, Eric W. "Ejemplo de cálculo de la varianza" . MathWorld .