Ver el código de «Algoritmos para calcular la varianza»

Los '''algoritmos para calcular la varianza''' juegan un papel importante en la [[estadística informática]]. Una dificultad clave en el diseño de buenos [[algoritmo]]s para este problema es que las fórmulas de la [[varianza]] pueden incluir sumas de cuadrados, lo que puede llevar a problemas de [[estabilidad numérica|inestabilidad numérica]] así como a [[rebosamiento]] cuando se trata de valores grandes.

== Algoritmo "ingenuo" ==
Una fórmula para calcular la varianza de una [[Población estadística|población]] completa de tamaño ''N'' es:

: <math>\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}. \!</math>

Usando la [[corrección de Bessel]] para calcular una estimación [[Sesgo estadístico|no sesgada]] de la varianza poblacional de una [[Muestra estadística|muestra]] finita de observaciones ''n'', la fórmula es:

: <math>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}. \!</math>

Por lo tanto, un algoritmo ''"ingenuo"'' para calcular la varianza estimada viene dado por el procedimiento:

<div style="margin-left: 35px; width: 600px">
{{framebox|blue}}
* Let {{math|''n'' ← 0, Sum ← 0, SumSq ← 0}}
* For each datum {{mvar|x}}:
** {{math|''n'' ← ''n'' + 1}}
** {{math|Sum ← Sum + ''x''}}
** {{math|SumSq ← SumSq + ''x'' × ''x''}}
* {{math|Var {{=}} (SumSq − (Sum × Sum) / n) / (n − 1)}}
{{frame-footer}}
</div>

Este algoritmo se puede adaptar fácilmente para calcular la varianza de una población finita: simplemente, basta con dividr entre ''N'' en lugar de entre ''n''&nbsp;-&nbsp;1 en la última línea.

Como {{math|SumSq}} y {{math|(Sum×Sum)/''n''}} pueden ser números muy similares, [[pérdida de significancia|cancelación]] de los restos puede llevar a que la [[Cifras significativas|precisión]] del resultado sea mucho menor que la precisión inherente de [[coma flotante]] utilizada para realizar el cálculo. Por lo tanto, se ha propuesto este algoritmo no debe utilizarse en la práctica.<ref name="Einarsson2005">{{cite book|first=Bo|last=Einarsson|title=Accuracy and Reliability in Scientific Computing|url=https://books.google.com/books?id=8hrDV5EbrEsC|accessdate=17 de febrero de 2013|date=1 de agosto de 2005|publisher=SIAM|isbn=978-0-89871-584-2|page=47}}</ref><ref name="Chan1983">{{cite journal|url=http://cpsc.yale.edu/sites/default/files/files/tr222.pdf|first1=Tony F.|last1=Chan|author1-link=Tony F. Chan|first2=Gene H.|last2=Golub|author2-link=Gene H. Golub|first3=Randall J.|last3=LeVeque|title=Algorithms for computing the sample variance: Analysis and recommendations|journal=The American Statistician|volume=37|issue=3|pages=242–247|year=1983|jstor=2683386}}</ref> y varios algoritmos alternativos, numéricamente estables.<ref name=":1">{{Cite journal|last=Schubert|first=Erich|last2=Gertz|first2=Michael|date=9 de julio de 2018|title=Numerically stable parallel computation of (co-)variance|url=http://dl.acm.org/citation.cfm?id=3221269.3223036|publisher=ACM|pages=10|doi=10.1145/3221269.3223036|isbn=9781450365055}}</ref> Esto es particularmente malo si la desviación estándar es pequeña en relación con la media. Sin embargo, el algoritmo puede mejorarse adoptando el método de la [[media falsa]].

=== Cálculo de datos desplazados ===

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

: <math>\operatorname{Var}(X-K)=\operatorname{Var}(X).</math>

Siendo <math>K</math> una constante cualquiera, se obtiene la nueva fórmula:

: <math>s^2 = \displaystyle\frac {\sum_{i=1}^n (x_i-K)^2 - (\sum_{i=1}^n (x_i-K))^2/n}{n-1}. \!</math>

Cuanto más cerca esté <math>K</math> del valor medio, más preciso será el resultado, pero simplemente eligiendo un valor dentro del rango de muestras se garantizará la estabilidad deseada. Si los valores <math>(x_i - K)</math> 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 en la fórmula es siempre más pequeño que el primero, por lo que no puede producirse ninguna cancelación.<ref name="Chan1983" />

Si se toma solo la primera muestra como <math>K</math>, el algoritmo se puede escribir en el [[Python|lenguaje de programación Python]] como

<syntaxhighlight lang="python">
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
</syntaxhighlight>
esta fórmula facilita también el cálculo incremental, que puede expresarse como
<syntaxhighlight lang="python">
K = n = Ex = Ex2 = 0.0

def add_variable (x):
     si (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 ():
     volver K + Ex / n

def get_variance ():
     retorno (Ex2 - (Ex * Ex) / n) / (n-1)
</syntaxhighlight>

== Algoritmo de dos pasos ==
Una aproximación alternativa, usando una fórmula diferente para la varianza, computa primero la media de la muestra,
:<math>\bar x = \frac {\sum_{j=1}^n x_j} n,</math>
y después computa la suma de los cuadrados de las diferencias con la media,
:<math>\text{sample variance} = s^2 = \displaystyle\frac {\sum_{i=1}^n (x_i - \bar x)^2}{n-1}, \!</math>
donde ''s'' es la desviación estándar. Este planteamiento viene dado por el siguiente [[pseudocódigo]]:
<syntaxhighlight lang="python">
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
</syntaxhighlight>

Este algoritmo es numéricamente estable si ''n'' es pequeño.<ref name="Einarsson2005"/><ref>{{cite book|first=Nicholas | last=Higham |title=Accuracy and Stability of Numerical Algorithms (2 ed) (Problem 1.10)| publisher=SIAM|year=2002}}</ref> Sin embargo, los resultados de estos dos algoritmos simples (''"ingenuo"'' y ''"de dos pasos"'') pueden depender excesivamente del orden de los datos y pueden dar resultados deficientes para conjuntos de datos muy grandes, debido a repetidos errores de redondeo en la acumulación de las sumas. Técnicas como la [[suma compensada]] se pueden usar para combatir este error hasta cierto punto.

== Algoritmo en línea de Welford ==
A menudo es útil poder calcular la variación en una sola pasada, inspeccionando cada valor <math>x_i</math> solo una vez; por ejemplo, cuando los datos se recopilan sin el almacenamiento suficiente para mantener todos los valores, o cuando los costos de acceso a la memoria dominan a los de cálculo. Para este [[algoritmo en línea]], se utiliza 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 [[media (matemáticas)|media]] y la varianza (estimada) de la secuencia, para un elemento adicional ''x''<sub>''n''</sub>. Aquí, {{overline|''x''}}<sub>''n''</sub> denota la media muestral de las primeras ''n'' muestras (''x''<sub>1</sub>, ..., ''x<sub>n</sub>''), ''s''{{supsub|2|''n''}} su varianza muestral, y ''σ''{{supsub|2|''n''}} su varianza poblacional.

:<math>\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} \!</math>

:<math>s^2_n = \frac{n-2}{n-1} \, s^2_{n-1} + \frac{(x_n - \bar x_{n-1})^2}{n} = s^2_{n-1} + \frac{(x_n - \bar x_{n-1})^2}{n} - \frac{s^2_{n-1}}{n-1}, \quad n>1 </math>

:<math>\sigma^2_n = \frac{(n-1) \, \sigma^2_{n-1} + (x_n - \bar x_{n-1})(x_n - \bar x_n)}{n} = \sigma^2_{n-1} + \frac{(x_n - \bar x_{n-1})(x_n - \bar x_n) - \sigma^2_{n-1}}{n}.</math>

Estas fórmulas sufren de inestabilidad numérica, ya que restan repetidamente un número pequeño de un número grande que se escala con ''n''. Una mejor cantidad para actualizar es la suma de cuadrados de las diferencias de la media actual, <math>\textstyle\sum_{i=1}^n (x_i - \bar x_n)^2</math>, que aquí denota <math>M_{2,n}</math>:

: <math>
\begin{align}
M_{2,n} & = M_{2,n-1} + (x_n - \bar x_{n-1})(x_n - \bar x_n) \\[4pt]
s^2_n & = \frac{M_{2,n}}{n} \\[4pt]
\sigma^2_n & = \frac{M_{2,n}}{n-1}
\end{align}
</math>

Este algoritmo fue ideado por Welford,<ref>{{cite journal |first=B. P. |last=Welford |year=1962 |title=Note on a method for calculating corrected sums of squares and products |journal=[[Technometrics]] |volume=4 |issue=3 |pages=419–420 |jstor=1266577 |doi=10.2307/1266577}}</ref><ref>[[Donald Knuth]] (1998). ''[[The Art of Computer Programming]]'', volume 2: ''Seminumerical Algorithms'', 3rd edn., p.&nbsp;232. Boston: Addison-Wesley.</ref> y se ha analizado exhaustivamente.<ref name="Chan1983" /><ref>{{cite journal |last=Ling |first=Robert F. |year=1974 |title=Comparison of Several Algorithms for Computing Sample Means and Variances |journal=Journal of the American Statistical Association |volume=69 |issue=348 |pages=859–866 |doi=10.2307/2286154}}</ref> También es común denotar <math>M_k = \bar x_k</math> y <math>S_k = M_{2,k}</math>.<ref>http://www.johndcook.com/standard_deviation.html</ref>

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

<syntaxhighlight lang="python">
# 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)
</syntaxhighlight>

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 incluida dentro del bucle. Para disponer de un algoritmo de dos pasos particularmente robusto con el fin de calcular la varianza, primero se puede calcular y restar una estimación de la media y luego usar este algoritmo con los residuos.

El ''algoritmo paralelo'' mostrado más adelante ilustra cómo combinar varios conjuntos de estadísticas calculadas en línea.

== Algoritmo incremental ponderado ==
El algoritmo puede extenderse para manejar pesos de muestra desiguales, reemplazando el simple contador ''n'' con la suma de pesos vista hasta ahora. West (1979)<ref>{{cite journal |first=D. H. D. |last=West |year=1979 |title=Updating Mean and Variance Estimates: An Improved Method |journal=[[Communications of the ACM]] |volume=22 |issue=9 |pages=532–535 |doi=10.1145/359146.359153}}</ref> sugiere este [[algoritmo incremental]]:

<syntaxhighlight lang="python">
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)
</syntaxhighlight>

== Algoritmo "paralelo" ==
Chan et al.<ref name=":0">{{Citation
  | last1 = Chan     | first1 = Tony F.      | author1-link = Tony F. Chan
  | last2 = Golub    | first2 = Gene H.      | author2-link = Gene H. Golub
  | last3 = LeVeque  | first3 = Randall J.
  | contribution = Updating Formulae and a Pairwise Algorithm for Computing Sample Variances.
  | title = Technical Report STAN-CS-79-773
  | publisher = Department of Computer Science, Stanford University
  | year = 1979
  | contribution-url =http://i.stanford.edu/pub/cstr/reports/cs/tr/79/773/CS-TR-79-773.pdf }}.</ref> observaron que el algoritmo en línea de Welford detallado anteriormente es un caso especial de un algoritmo que funciona para cualquier partición de la muestra <math>X</math> en los conjuntos <math>X_A</math>, <math>X_B</math>:
: <math>\delta\! = \bar x_B - \bar x_A</math>
: <math>\bar x_X = \bar x_A + \delta\cdot\frac{n_B}{n_X}</math>
: <math>M_{2,X} = M_{2,A} + M_{2,B} + \delta^2\cdot\frac{n_A n_B}{n_X}</math>.
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 <math>n_A \approx n_B</math> y ambos son grandes, porque el error numérico en <math>\bar x_B - \bar x_A</math> no se reduce de la forma en que se encuentra en el caso <math>n_B = 1</math>. En tales casos, es preferible <math>\bar x_X = \frac{n_A \bar x_A + n_B \bar x_B}{n_A + n_B}</math>.<syntaxhighlight lang="python">
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)
</syntaxhighlight> Esto se puede generalizar para permitir la ''paralelización'' con ''[[Extensiones Vectoriales Avanzadas]]'' (EVA), con [[Unidad de procesamiento gráfico|unidades de procesamiento gráfico]] y [[Clúster (informática)|agregados de computadoras]], y para la covarianza.<ref name=":1" />

== Ejemplo ==
Supóngase que todas las operaciones de punto flotante utilizan aritmética [[IEEE coma flotante|IEEE 754 de doble precisión]] estándar. Considérese la muestra (4, 7, 13, 16) de una población infinita. En base a 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 pasos calculan estos valores correctamente.

A continuación, considérese la muestra ({{nowrap|10<sup>8</sup>&nbsp;+&nbsp;4}}, {{nowrap|10<sup>8</sup>&nbsp;+&nbsp;7}}, {{nowrap|10<sup>8</sup>&nbsp;+&nbsp;13}}, {{nowrap|10<sup>8</sup>&nbsp;+&nbsp;16}}), 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.33333333333333332 en lugar de 30.

Si bien esta pérdida de precisión puede ser tolerable y vista como un defecto menor del algoritmo "ingenuo", un aumento adicional del desplazamiento hace que el error sea catastrófico. Considérese la muestra ({{nowrap|10<sup>9</sup>&nbsp;+&nbsp;4}}, {{nowrap|10<sup>9</sup>&nbsp;+&nbsp;7}}, {{nowrap|10<sup>9</sup>&nbsp;+&nbsp;13}}, {{nowrap|10<sup>9</sup>&nbsp;+&nbsp;16}}). De nuevo, la varianza estimada de la población de 30 se calcula correctamente mediante el algoritmo de dos pasos, pero el algoritmo "ingenuo" ahora lo calcula como −170.66666666666666. Este es un problema serio 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<ref>{{Citation
 | last=Terriberry
 | first=Timothy B.
 | year=2007
 | title=Computing Higher-Order Moments Online
 | url=http://people.xiph.org/~tterribe/notes/homs.html
 | fechaacceso=9 de abril de 2019
 | urlarchivo=https://web.archive.org/web/20140423031833/http://people.xiph.org/~tterribe/notes/homs.html
 | fechaarchivo=23 de abril de 2014
 }}</ref> amplía las fórmulas de Chan para calcular el tercer y cuarto [[momento central]], necesarios, por ejemplo, al estimar la [[asimetría estadística]] y la [[curtosis]]:
: <math>
\begin{align}
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_AM_{2,B} - n_BM_{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_AM_{3,B} - n_BM_{3,A}}{n_X}
\end{align}</math>

Aquí, los <math>M_k</math> son nuevamente las sumas de las potencias de las diferencias conrespecto a la media <math>\sum(x - \overline{x})^k</math>, dando
: <math>
\begin{align}
& \text{skewness} = g_1 = \frac{\sqrt{n} M_3}{M_2^{3/2}}, \\[4pt]
& \text{kurtosis} = g_2 = \frac{n M_4}{M_2^2}-3.
\end{align}
</math>

Para el caso incremental (es decir, <math>B = \{x\}</math>), esto se simplifica a:
: <math>
\begin{align}
\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{align}
</math>

Al preservar el valor <math>\delta / n</math>, solo se necesita una operación de división y, por lo tanto, las estadísticas de orden superior se pueden calcular con un reducido costo incremental.

Un ejemplo del algoritmo en línea para la curtosis implementado como se describe es:
<syntaxhighlight lang="python">
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
</syntaxhighlight>

Pébaÿ<ref>{{Citation
 | last=Pébaÿ
 | first=Philippe
 | year=2008
 | contribution=Formulas for Robust, One-Pass Parallel Computation of Covariances and Arbitrary-Order Statistical Moments
 | title=Technical Report SAND2008-6212
 | publisher=Sandia National Laboratories
 | contribution-url=http://infoserve.sandia.gov/sand_doc/2008/086212.pdf
 }}</ref>
extiende aún más estos resultados a [[momento central|momentos centrales]] de orden arbitrario, para los casos incrementales y por pares, y posteriormente Pébaÿ et al.<ref>{{Citation
 | last1=Pébaÿ
 | first1=Philippe
 | last2=Terriberry
 | first2=Timothy
 | last3=Kolla
 | first3=Hemanth
 | last4=Bennett
 | first4=Janine
 | year=2016
 | title=Numerically Stable, Scalable Formulas for Parallel and Online Computation of Higher-Order Multivariate Central Moments with Arbitrary Weights
 | journal=Computational Statistics
 | volume=31
 | issue=4
 | pages=1305–1325
 | publisher=Springer
 | doi=10.1007/s00180-015-0637-z
 }}</ref>
lo hicieron para momentos ponderados y compuestos. También se pueden encontrar fórmulas similares para la [[covarianza]].

Choi y Sweetman<ref name="Choi2010">{{Citation
 |last1  = Choi
 |first1 = Myoungkeun
 |last2  = Sweetman
 |first2 = Bert
 |year   = 2010
 |title  = Efficient Calculation of Statistical Moments for Structural Health Monitoring
 |url    = https://www.semanticscholar.org/paper/Efficient-Calculation-of-Statistical-Moments-for-Choi-Sweetman/91379baabdc47568327f3646a587af5ee5bb7698
 |journal= Journal of Structural Health Monitoring
 |volume=9
 |issue=1
 |pages=13–24
 |doi=10.1177/1475921709341014
}}</ref>
ofrecieron dos métodos alternativos para calcular la asimetría y la curtosis, cada uno de los cuales puede ahorrar importantes requisitos 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, que efectivamente se convierte en un [[algoritmo de un paso]] para momentos más altos. Un beneficio es que los cálculos estadísticos del momento se pueden realizar con una [[precisión arbitraria]], de manera que se pueden 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 manera convencional: el rango de valores potenciales es dividido en contenedores y la cantidad de ocurrencias dentro de cada contenedor se recuenta y se traza de manera tal que el área de cada rectángulo es igual a la porción de los valores de muestra dentro de ese contenedor:

: <math> H(x_k)=\frac{h(x_k)}{A}</math>

donde <math>h(x_k)</math> y <math>H(x_k)</math> representan la frecuencia y la frecuencia relativa en el intervalo <math>x_k</math> y <math>A= \sum_{k=1}^K h(x_k) \,\Delta x_k</math> es el área total del histograma. Después de esta normalización, los momentos en bruto <math>n</math> y los momentos centrales de <math>x(t)</math> se pueden calcular a partir del histograma relativo:

: <math>
 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
</math>

: <math>
 \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
</math>

donde el superíndice <math>^{(h)}</math> indica que los momentos se calculan a partir del histograma. Para un ancho constante del contenedor <math>\Delta x_k=\Delta x</math> estas dos expresiones se pueden simplificar utilizando <math>I= A/\Delta x</math>:

: <math>
 m_n^{(h)}= \frac{1}{I} \sum_{k=1}^K x_k^n \, h(x_k)
</math>

: <math>
 \theta_n^{(h)}= \frac{1}{I} \sum_{k=1}^K \Big(x_k-m_1^{(h)}\Big)^n h(x_k)
</math>

El segundo enfoque de Choi y Sweetman<ref name="Choi2010" /> es una metodología analítica para combinar momentos estadísticos de segmentos individuales de un historial en el tiempo, de modo que los momentos globales resultantes sean los del historial del tiempo completo. 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 computados en tiempos secuenciales.

Si se conocen conjuntos <math>Q</math> de momentos estadísticos:
<math>(\gamma_{0,q},\mu_{q},\sigma^2_{q},\alpha_{3,q},\alpha_{4,q})
\quad </math> para <math>q=1,2,\ldots,Q </math>, entonces cada <math>\gamma_n</math> puede
expresarse en términos de los momentos en bruto <math>n</math> equivalentes:

: <math>
\gamma_{n,q}= m_{n,q} \gamma_{0,q} \qquad \quad \textrm{para } \quad n=1,2,3,4  \quad \text{ and } \quad q = 1,2, \dots ,Q
</math>

donde <math>\gamma_{0,q}</math> generalmente se toma como la duración del historial de tiempo de <math>q^{th}</math>, o el número de puntos si <math>\Delta t</math> es constante.

La ventaja de expresar los momentos estadísticos en términos de <math>\gamma</math> es que los conjuntos de <math>Q</math> se pueden combinar por adición, y no hay límite superior en el valor de <math>Q</math>.

: <math>
 \gamma_{n,c}= \sum_{q=1}^Q \gamma_{n,q} \quad \quad \text{para } n=0,1,2,3,4
</math>

donde el subíndice <math>_c</math> representa el historial en el tiempo concatenado o <math>\gamma</math> combinado. Estos valores combinados de <math>\gamma</math> se pueden transformar inversamente en momentos sin procesar que representan el historial en el tiempo concatenado completo

: <math>
 m_{n,c}=\frac{\gamma_{n,c}}{\gamma_{0,c}} \quad \text{para } n=1,2,3,4
</math>

Relaciones conocidas entre los momentos sin procesar (<math>m_n</math>) y los momentos centrales (<math> \theta_n = \operatorname E[(x-\mu)^n])</math>)
se utilizan después para calcular los momentos centrales de la historia temporal concatenada. Finalmente, los momentos estadísticos de la historia concatenada se computan a partir de los momentos centrales:

: <math>
 \mu_c=m_{1,c}
 \qquad \sigma^2_c=\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
</math>

== Covarianza ==
Se pueden usar algoritmos muy similares para calcular el [[covarianza]].

=== Algoritmo "ingenuo" ===
El algoritmo "ingenuo" es:
: <math>\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}. \!</math>

Para el algoritmo anterior, se podría usar el siguiente código de Python:
<syntaxhighlight lang="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
</syntaxhighlight>

=== Con estimación de la media ===
En cuanto a la varianza, la covarianza de dos variables aleatorias también es invariante al cambio, por lo que dados dos valores constantes <math>k_x</math> y <math>k_y,</math> pueden escribirse:

: <math>\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}. \!</math>

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

<syntaxhighlight lang="python">
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
</syntaxhighlight>

=== Algoritmo de dos pasos ===
El algoritmo de dos pasos primero calcula las medias muestrales y luego la covarianza:
: <math>\bar x = \displaystyle \sum_{i=1}^n x_i/n</math>
: <math>\bar y = \displaystyle \sum_{i=1}^n y_i/n</math>
: <math>\operatorname{Cov}(X,Y) = \displaystyle\frac {\sum_{i=1}^n (x_i - \bar x)(y_i - \bar y)}{n}. \!</math>

Se puede escribir como:
<syntaxhighlight lang="python">
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
</syntaxhighlight>

Una versión compensada ligeramente más precisa realiza el algoritmo "ingenuo" completamente en los residuos. Las sumas finales <math>\textstyle\sum x_i</math> y <math>\textstyle\sum y_i</math> '' deberían ''ser cero'', pero la segunda pasada compensa cualquier pequeño error.

=== Algoritmo en línea ===
Existe un algoritmo estable de una pasada, similar al algoritmo en línea para calcular la varianza, que calcula el momento <math>\textstyle C_n = \sum_{i=1}^n (x_i - \bar x_n)(y_i - \bar y_n)</math>:
: <math>\begin{alignat}{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{alignat}</math>
La asimetría aparente en esa última ecuación se debe al hecho de que <math>\textstyle (x_n - \bar x_n) = \frac{n-1}{n}(x_n - \bar x_{n-1})</math>, por lo que ambos términos de actualización son iguales a <math>\textstyle \frac{n-1}{n}(x_n - \bar x_{n-1})(y_n - \bar y_{n-1})</math>. Se puede lograr una mayor precisión si primero se calculan las medias y luego se usa el algoritmo estable de una pasada en los residuos.

Así se puede calcular la covarianza como
: <math>\begin{align}
\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) + (y_n - \bar y_n)(x_n - \bar x_{n-1})}{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}.
\end{align}</math>

<syntaxhighlight lang="python">
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)
</syntaxhighlight>

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

<syntaxhighlight lang="python">
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)
</syntaxhighlight>

Del mismo modo, existe una fórmula para combinar las covarianzas de dos conjuntos que se pueden usar para paralelizar el cálculo:<ref name=":1" />

: <math>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}. </math>

=== Versión ponderada por lotes ===
También existe una versión del algoritmo en línea ponderado que se actualiza por lotes: haciendo que
<math>w_1, \dots w_N</math> denote los pesos, se puede escribir

: <math>\begin{alignat}{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{alignat}</math>

La covarianza se puede calcular como

: <math>\operatorname{Cov}_N(X,Y) = \frac{C_N}{\sum_{i=1}^{N} w_i}</math>

== Véase también ==
* [[Fórmula algebraica de la varianza]]
* [[Algoritmo sumatorio de Kahan]]
* [[Desviaciones cuadráticas respecto a la media]]

== Referencias ==
{{listaref}}

== Enlaces externos ==
* {{MathWorld|title=Sample Variance Computation|urlname=SampleVarianceComputation}}

{{Control de autoridades}}
[[Categoría:Dispersión estadística]]
[[Categoría:Artículos con pseudocódigos de ejemplo]]
[[Categoría:Algoritmos estadísticos]]