Los algoritmos para el cálculo de la varianza - Algorithms for calculating variance


De Wikipedia, la enciclopedia libre

Los algoritmos para el cálculo de la varianza juegan un papel importante en la estadística computacional . Una de las principales dificultades en el diseño de buenos algoritmos para este problema es que las fórmulas de la varianza pueden implicar sumas de cuadrados, que pueden conducir a la inestabilidad numérica , así como al desbordamiento aritmético cuando se trata de grandes valores.

algoritmo ingenuo

Una fórmula para el cálculo de la varianza de una entera población de tamaño N es:

El uso de la corrección de Bessel para calcular un imparcial estimación de la varianza de la población a partir de un número finito de muestra de n observaciones, la fórmula es:

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

  • Deje n ← 0, 0 ← Suma, SUMSQ ← 0
  • Para cada dato x :
    • nn + 1
    • ← Suma Suma + x
    • SUMSQ ← SUMSQ + x x x
  • Var = (SUMSQ - (Sum × Sum) / n) / (n - 1)

Este algoritmo se puede adaptar fácilmente para calcular la varianza de una población finita: simplemente dividir por N en lugar de n  - 1 en la última línea.

Debido SUMSQ y (Sum Sum ×) / n pueden ser números muy similares, la cancelación puede conducir a la precisión del resultado ser mucho menor que la precisión inherente de la aritmética de punto flotante utilizada para realizar el cálculo. Así, este algoritmo no debe utilizarse en la práctica, y varios alternativo, estable numéricamente, se han propuesto algoritmos. 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 mediante la adopción del método de la media asumida .

Computación cambió de datos

Podemos utilizar una propiedad de la varianza para evitar la cancelación catastrófica en esta fórmula, es decir, la varianza es invariante con respecto a los cambios en un parámetro de localización

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

cuanto más cerca es el valor medio más preciso será el resultado será, pero sólo la elección de un valor dentro de las muestras varían garantice la estabilidad deseada. Si los valores son pequeños entonces no hay ningún problema 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 por lo tanto se puede producir el primero no cancelación.

Si tomamos sólo la primera muestra que el algoritmo puede ser escrito en el lenguaje de programación Python como

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 facilita así el cálculo incrementales, que puede ser expresado como

K = n = Ex = Ex2 = 0.0

def add_variable(x):
    if (n == 0):
      K = x
    n = n + 1
    Ex += x - K
    Ex2 += (x - K) * (x - K)

def remove_variable(x):
    n = n - 1
    Ex -= (x - K)
    Ex2 -= (x - K) * (x - K)

def get_meanvalue():
    return K + Ex / n

def get_variance():
    return (Ex2 - (Ex*Ex)/n) / (n-1)

algoritmo de dos pasos

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

,

y después se calcula la suma de los cuadrados de las diferencias con respecto a la media,

,

donde s es la desviación estándar. Esto se da por la siguiente pseudocódigo:

def two_pass_variance(data):
    n = sum1 = sum2 = 0

    for x in data:
        n += 1
        sum1 += x

    mean = sum1 / n

    for x in data:
        sum2 += (x - mean)*(x - mean)

    variance = sum2 / (n - 1)
    return variance

Este algoritmo es numéricamente estable si n es pequeño. Sin embargo, los resultados de estos dos simples algoritmos ( "vírgenes" y "Two-pass") pueden depender excesivamente en el orden de los datos y pueden dar resultados pobres para grandes conjuntos de datos debido a un error de redondeo se repite en la acumulación de la sumas. Técnicas tales como la suma compensada se pueden utilizar para combatir este error en un grado.

algoritmo de línea de Welford

A menudo es útil ser capaz de calcular la varianza de una sola pasada, la inspección de cada valor sólo una vez; por ejemplo, cuando los datos están siendo recogidos y sin suficiente capacidad de almacenamiento para guardar todos los valores, o cuando los costos de acceso a la memoria dominar a los de computación. Para un tal algoritmo en línea , una relación de recurrencia se requiere entre las cantidades de las que las estadísticas requeridas se pueden calcular de una manera estable numéricamente.

Las siguientes fórmulas se pueden utilizar para actualizar la media y (estimado) varianza de la secuencia, para un adicional elemento x n . Aquí, x n denota la media de la muestra de las primeras n muestras ( x 1 , ..., x n ), es 2 n la varianza de la muestra, y σ 2 n la varianza de la población.

Estas fórmulas sufren de inestabilidad numérica, ya que restan en repetidas ocasiones un número pequeño de un gran número que escala con n . Una mejor cantidad para la actualización es la suma de los cuadrados de las diferencias respecto a la media actual, aquí denota :

Este algoritmo se encontró por Welford, y se ha analizado a fondo. También es común para denotar y .

Un ejemplo de implementación de Python para el algoritmo de Welford es la siguiente.

# for a new value newValue, compute the new count, new mean, the new M2.
# mean accumulates the mean of the entire dataset
# M2 aggregates the squared distance from the mean
# count aggregates the number of samples seen so far
def update(existingAggregate, newValue):
    (count, mean, M2) = existingAggregate
    count += 1 
    delta = newValue - mean
    mean += delta / count
    delta2 = newValue - mean
    M2 += delta * delta2

    return (count, mean, M2)

# retrieve the mean, variance and sample variance from an aggregate
def finalize(existingAggregate):
    (count, mean, M2) = existingAggregate
    (mean, variance, sampleVariance) = (mean, M2/count, M2/(count - 1)) 
    if count < 2:
        return float('nan')
    else:
        return (mean, variance, sampleVariance)

Este algoritmo es mucho menos propenso a la pérdida de precisión debido a la cancelación catastrófica , pero podría no ser tan eficiente debido a la operación de división dentro del bucle. Para una particularmente robusto algoritmo de dos pasos para el cálculo de la varianza, se puede primero 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 incremento ponderado

El algoritmo se puede ampliar para manejar pesos de muestra desiguales, en sustitución de la simple contador n con la suma de pesos visto hasta ahora. West (1979) sugiere que este algoritmo incrementales :

def weighted_incremental_variance(dataWeightPairs):
    wSum = wSum2 = mean = S = 0

    for x, w in dataWeightPairs:  # Alternatively "for x, w in zip(data, weights):"
        wSum = wSum + w
        wSum2 = wSum2 + w*w
        meanOld = mean
        mean = meanOld + (w / wSum) * (x - meanOld)
        S = S + w * (x - meanOld) * (x - mean)

    population_variance = S / wSum
    # Bessel's correction for weighted samples
    # Frequency weights
    sample_frequency_variance = S / (wSum - 1)
    # Reliability weights
    sample_reliability_variance = S / (wSum - wSum2/wSum)

algoritmo paralelo

Chan et al. en cuenta que el anterior algoritmo de "on-line" es un caso especial de un algoritmo que funciona para cualquier partición de la muestra en conjuntos , :

.

Esto puede ser útil cuando, por ejemplo, múltiples unidades de procesamiento pueden ser asignados a partes discretas de la entrada.

El método de Chan para la estimación de la media es numéricamente inestable cuando y ambos son grandes, debido a que el error numérico en el que no se escala en la forma en que lo es en el caso. En tales casos, prefiera .

def parallel_variance(avg_a, count_a, var_a, avg_b, count_b, var_b):
    delta = avg_b - avg_a
    m_a = var_a * (count_a - 1)
    m_b = var_b * (count_b - 1)
    M2 = m_a + m_b + delta ** 2 * count_a * count_b / (count_a + count_b)
    return M2 / (count_a + count_b - 1)

Ejemplo

Suponga que todas las operaciones de punto flotante uso estándar IEEE 754 de doble precisión aritmética. Considere la muestra (4, 7, 13, 16) de una población infinita. Sobre la base de esta muestra, la media de la población estimada es de 10, y la estimación no sesgada de la varianza de la población es de 30. Tanto el algoritmo ingenuo y de dos pasadas algoritmo de calcular estos valores correctamente.

Consideremos a continuación la muestra ( 10 8  + 4 , 10 8  + 7 , 10 8  + 13 , 10 8  + 16 ), que da lugar a la misma varianza estimada como la primera muestra. El algoritmo de dos pasos calcula esta estimación de la varianza correcta, pero el algoritmo ingenuo vuelve 29,333333333333332 en lugar de 30.

Si bien esta pérdida de precisión puede ser tolerable y visto como un defecto de menor importancia del algoritmo ingenuo, aumentando aún más el desplazamiento hace que el error catastrófico. Considere la muestra ( 10 9  + 4 , 10 9  + 7 , 10 9  + 13 , 10 9  + 16 ). De nuevo, la varianza de la población estimada de 30 se calcula correctamente por el algoritmo de dos pasos, pero el algoritmo ingenuo ahora calcula como -170.66666666666666. Este es un serio problema con el algoritmo ingenuo y se debe a la cancelación catastrófica en la resta de dos números similares en la etapa final del algoritmo.

estadísticas de orden superior

Terriberry extiende fórmulas de Chan para calcular el tercer y cuarto momentos centrales , necesarias, por ejemplo, al estimar la asimetría y curtosis :

Aquí el son de nuevo las sumas de las potencias de las diferencias respecto a la media , dando

oblicuidad:
curtosis:

Para el caso incremental (es decir, ), esto se simplifica a:

Al preservar el valor , sólo se necesita una operación de división y las estadísticas de orden superior por lo tanto puede ser calculada por poco coste incremental.

Un ejemplo del algoritmo de línea para kurtosis implementado como se describe es:

def online_kurtosis(data):
    n = mean = M2 = M3 = M4 = 0

    for x in data:
        n1 = n
        n = n + 1
        delta = x - mean
        delta_n = delta / n
        delta_n2 = delta_n * delta_n
        term1 = delta * delta_n * n1
        mean = mean + delta_n
        M4 = M4 + term1 * delta_n2 * (n*n - 3*n + 3) + 6 * delta_n2 * M2 - 4 * delta_n * M3
        M3 = M3 + term1 * delta_n * (n - 2) - 3 * delta_n * M2
        M2 = M2 + term1

    kurtosis = (n*M4) / (M2*M2) - 3
    return kurtosis

Pebay extiende aún más estos resultados a-orden arbitrario momentos centrales , para la incrementales y los casos por pares, y posteriormente Pebay et al. para momentos ponderados y compuestas. También se puede encontrar que hay fórmulas similares para covarianza .

Choi y Sweetman ofrecen dos métodos alternativos para calcular la asimetría y curtosis, cada uno de los cuales pueden ahorrar los requisitos de memoria de ordenador sustancial y tiempo de CPU en ciertas aplicaciones. El primer enfoque consiste en calcular los momentos estadísticos mediante la separación de los datos en los contenedores y luego calcular los momentos de la geometría de la histograma resultante, que se convierte efectivamente en un algoritmo de una sola pasada para momentos más elevados. Uno de los beneficios es que los cálculos estadísticos momento pueden llevarse a cabo con una precisión arbitraria de tal manera que los cálculos se pueden sintonizar a la precisión de, por ejemplo, el formato de almacenamiento de datos o el hardware de medición original. Un histograma relativa de una variable aleatoria puede ser construido de la manera convencional: el rango de valores posibles se divide en compartimientos y el número de ocurrencias dentro de cada bin se cuentan y se representa tal que el área de cada rectángulo es igual a la parte de los valores de la muestra dentro de ese bin:

donde y representar la frecuencia y la frecuencia relativa en bin y es el área total del histograma. Después de esta normalización, los momentos primas y momentos centrales de se pueden calcular a partir del histograma relativo:

donde el superíndice indica los momentos se calculan a partir del histograma. Para anchura bin constante estas dos expresiones se pueden simplificar usando :

El segundo enfoque de Choi y Sweetman es una metodología analítica combinar momentos estadísticos de segmentos individuales de un tiempo de la historia de tal manera que los momentos generales resultantes son los de la completa de historia de tiempo. Esta metodología podría ser utilizado 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 computados a veces secuenciales.

Si se conocen conjuntos de momentos estadísticos: para , a continuación, cada uno se puede expresar en términos de los equivalentes momentos primas:

donde se toma generalmente como la duración del tiempo de la historia, o el número de puntos si es constante.

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

donde el subíndice representa el tiempo-historia concatenada o combinados . Estos valores combinados de lata entonces ser transformados inversamente en momentos primas que representan la concatenado tiempo-historia completa

Relaciones conocidas entre los momentos crudos ( ) y los momentos centrales ( ) se utilizan para calcular los momentos centrales de la concatenados tiempo-historia. Por último, los momentos estadísticos de la historia concatenada se calculan a partir de los momentos centrales:

covarianza

Algoritmos muy similares se pueden utilizar para calcular la covarianza .

algoritmo ingenuo

El algoritmo es ingenuo:

Para el algoritmo anterior, se podría utilizar el siguiente código 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 se desplace invariante, por lo Dados dos valores constantes y se puede escribir:

y otra vez la elección de un valor dentro del rango de valores se estabilizará la fórmula contra la cancelación catastrófica, así como hacerlo más robusto frente a grandes sumas de dinero. Tomando el primer valor de cada conjunto de datos, el algoritmo se puede escribir como:

def shifted_data_covariance(dataX, dataY):
   n = len(dataX)
   if (n < 2):
     return 0
   kx = dataX[0]
   ky = dataY[0]
   Ex = Ey = Exy = 0
   for ix, iy in zip(dataX, dataY):
      Ex += ix - kx
      Ey += iy - ky
      Exy += (ix - kx) * (iy - ky)
   return (Exy - Ex * Ey / n) / n

De dos pasadas

El algoritmo de dos pasos calcula primero los medios de muestra, y luego la covarianza:

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 ligeramente más precisa compensado realiza el algoritmo ingenuo completa de los residuos. Las sumas finales y deben ser cero, pero el segundo pase compensa cualquier pequeño error.

En línea

Existe un algoritmo de un paso estable, similar al algoritmo en línea para el cálculo de la varianza, que calcula co-momento :

La asimetría de manifiesto en la última ecuación es debido al hecho de que , por lo tanto en términos de actualización son iguales a . Incluso una mayor precisión se puede lograr mediante el cálculo de los medios primero, a continuación, utilizando el algoritmo de un paso estable en los residuos.

Por lo tanto podemos calcular la covarianza como

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 podemos hacer 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)

Del mismo modo, existe una fórmula para combinar las covarianzas de dos conjuntos que se pueden utilizar para paralelizar el cálculo:

versión ponderada por lotes

Una versión del algoritmo en línea ponderada que no actualiza lotes cor también existe: dejar que denotan los pesos, podemos escribir

La covarianza entonces puede calcularse como

Ver también

referencias

  1. ^ Un b Bo Einarsson (1 de agosto de 2005). La precisión y la fiabilidad de Computación Científica . SIAM. pag. 47. ISBN  978-0-89871-584-2 . Consultado el 17 de de febrero de 2013 .
  2. ^ Un b T.F.Chan, GH Golub y RJ LeVeque (1983). " " Algoritmos para el cálculo de la varianza de la muestra: Análisis y recomendaciones", el estadístico de América, 37" (PDF) : 242-247.
  3. ^ Un b Schubert, Erich; Gertz, Michael (09/07/2018). "Computación paralela numéricamente estable de (co) varianza" . ACM: 10. doi : 10.1145 / 3.221.269,3223036 . ISBN  9781450365055 .
  4. ^ Higham, Nicholas (2002). Precisión y estabilidad de Algoritmos Numéricos (2 ed) (Problema 1.10) . SIAM.
  5. ^ BP Welford (1962). "Nota sobre un método para el cálculo de sumas corregidas de cuadrados y productos" . Technometrics 4 (3): 419-420.
  6. ^ Donald E. Knuth (1998). The Art of Computer Programming , volumen 2: Seminumerical Algoritmos ., 3ª edición, p. 232. Boston: Addison-Wesley.
  7. ^ Chan, de Tony F .; Golub, Gene H .; LeVeque, Randall J. (1983). Algoritmos para calcular la varianza de la muestra: Análisis y Recomendaciones. El Estadístico de América 37, 242-247. https://www.jstor.org/stable/2683386
  8. ^ Ling, Robert F. (1974). Comparación de varios algoritmos para medios de cálculo de ejemplo y varianzas. Revista de la Asociación Americana de Estadística, vol. 69, No. 348, 859-866. doi : 10.2307 / 2286154
  9. ^ http://www.johndcook.com/standard_deviation.html
  10. ^ DHD West (1979). Comunicaciones de la ACM , 22, 9, 532-535: actualización de la estimación media y la varianza: Un método mejorado
  11. ^ Chan, de Tony F. ; Golub, Gene H. ; LeVeque, Randall J. (1979), "Fórmulas Actualizando y un algoritmo de Pairwise para calcular las varianzas de las muestras." (PDF) , Informe Técnico STAN-CS-79-773 , del Departamento de Ciencias de la Computación de la Universidad de Stanford.
  12. ^ Terriberry, Timothy B. (2007), Informática de orden superior momentos en línea
  13. ^ Pebay, Philippe (2008), "Fórmulas para robusta de una pasada Computación, en paralelo de covarianzas y momentos estadísticos arbitrarias-Order" (PDF) , Informe Técnico SAND2008-6212 , Sandia National Laboratories
  14. ^ Pebay, Philippe; Terriberry, Timothy; Kolla, Hemanth; Bennett, Janine (2016), "Fórmulas escalables, estable numéricamente para paralelo y en línea Cálculo de orden superior Momentos multivariado centrales con pesos arbitraria" , Estadística Computacional , Springer
  15. ^ Un b Choi, Myoungkeun; Sweetman, Bert (2010), "Cálculo eficiente de los momentos estadísticos para supervisar la salud estructural" , Journal of Structural Health Monitoring , 9 (1)

enlaces externos