Algoritmos para el cálculo de la varianza

Los algoritmos para calcular la varianza desempeñan un papel fundamental en la estadística computacional . Una dificultad clave en el diseño de buenos algoritmos para este problema radica en que las fórmulas para la varianza pueden incluir sumas de cuadrados, lo que puede provocar inestabilidad numérica y desbordamiento aritmético al trabajar con valores grandes.

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

${\displaystyle \sigma ^{2}={\overline {(x-{\bar {x}})^{2}}}={\overline {(x^{2})}}-{\bar {x}}^{2}={\frac {\sum _{i=1}^{N}x_{i}^{2}}{N}}-\left({\frac {\sum _{i=1}^{N}x_{i}}{N}}\right)^{2}}$

Utilizando la corrección de Bessel para calcular una estimación insesgada de la varianza poblacional 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 provocar 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 problemático 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 usarse para evitar la cancelación catastrófica en esta fórmula.

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

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

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

cuanto más cerca${\displaystyle K}$cuanto más cercano al valor medio, más preciso será el resultado, pero simplemente elegir un valor dentro del rango de las muestras garantizará la estabilidad deseada. Si los valores${\displaystyle (x_{i}-K)}$Si 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 siempre es menor que el primero, por lo que no puede haber cancelación. ^{[ 2 ]}

Si solo se toma la primera muestra como${\displaystyle K}$El algoritmo se puede escribir 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 += 1 Ex += x - K Ex2 += ( x - K ) ** 2 varianza = ( Ex2 - Ex ** 2 / n ) / ( n - 1 ) # usa n en lugar de (n-1) si quieres calcular la varianza exacta de los datos dados # usa (n-1) si los datos son muestras de una población mayor return varianza

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

A menudo resulta útil poder calcular la varianza en una sola pasada , inspeccionando cada valor.${\displaystyle x_{i}}$Solo una vez; por ejemplo, cuando se recopilan datos sin suficiente espacio de almacenamiento para guardar todos los valores, o cuando los costos de acceso a la memoria superan a los de computación. Para un algoritmo en línea de este tipo , se requiere una relación de recurrencia entre las cantidades a partir de la cual se puedan calcular las estadísticas necesarias de forma numéricamente estable.

clase ShiftDataVariance : def __init __ ( self ) : self.K = 0.0 self.n = 0 self.Ex = 0.0 self.Ex2 = 0.0​​​​​​def add_variable ( self , x : float ) : if self.n == 0 : self.K = x self.n + = 1 self.Ex + = x - self.K self.Ex2 + = ( x - self.K ) ** 2​​​​​​​​​​def remove_variable ( self , x : float ) : self.n - = 1 self.Ex - = x - self.K self.Ex2 - = ( x - self.K ) ** 2​​​​​​def add_variables ( self , xs : list [ float ] ) : # Python usa el algoritmo de Neumaier para la función integrada sum(), # que es más preciso que un simple bucle . if self.n == 0 and xs : self.K = xs [ 0 ] self.n + = len ( xs ) self.Ex + = sum ( x - self.K for x in xs ) self.Ex2 + = sum ( ( x - self.K ) ** 2 for x in xs )def obtener_media ( self ) -> float : return self . K + self . Ex / self . ndef get_variance ( self ) -> float : return ( self . Ex2 - self . Ex ** 2 / self . n ) / ( self . n - 1 )

clase WelfordVariance : def __init__ ( self ): # Comparación con ShiftDataVariance: self.mean = 0.0 # = K + Ex / n self.count = 0 # = n self.M2 = 0.0 # = Ex2 - (Ex)^2 / ndef add_variable ( self , x : float ): self . count += 1 old_mean = self . mean self . mean += ( x - self . mean ) / self . count self . M2 += ( x - old_mean ) * ( x - self . mean )def remove_variable ( self , x : float ): self . count -= 1 new_mean = self . mean self . mean -= ( x - self . mean ) / self . count self . M2 -= ( x - new_mean ) * ( x - self . mean )def obtener_media ( self ) -> float : return self . mediadef obtener_varianza ( self ) -> float : return self . M2 / self . countdef get_sample_variance ( self ) -> float : return self . M2 / ( self . count - 1 )

from collections import namedtuple WeightedVariances = namedtuple ( "WeightedVariances" , "pop freq reli" ) def weighted_incremental_variance ( data_weight_pairs ): w_sum = w_sum2 = mean = S = 0para x , w en data_weight_pairs : w_sum = w_sum + w w_sum2 = w_sum2 + w ** 2 media_antigua = media media = media_antigua + ( w / w_sum ) * ( x - media_antigua ) S = S + w * ( x - media_antigua ) * ( x - media )varianza_poblacional = S / w_sum # Corrección de Bessel para muestras ponderadas # Ponderaciones de frecuencia varianza_frecuencia_muestra = S / ( w_sum - 1 )# Pesos de confiabilidad varianza_confiabilidad_muestra = S / ( w_sum - w_sum2 / w_sum ) return VarianzasPonderadas ( varianza_población , varianza_frecuencia_muestra , varianza_confiabilidad_muestra )

Chan et al. ^{[ 11 ]} 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.${\displaystyle A}$y${\displaystyle 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 resultar útil cuando, por ejemplo, se pueden asignar varias unidades de procesamiento a partes discretas de la entrada.

El método de Chan para estimar la media es numéricamente inestable cuando${\displaystyle n_{A}\approx n_{B}}$y ambos son grandes, debido al error numérico en${\displaystyle \delta ={\bar {x}}_{B}-{\bar {x}}_{A}}$no se reduce de la misma manera que en el${\displaystyle n_{B}=1}$caso. En tales casos, prefiera${\textstyle {\bar {x}}_{AB}={\frac {n_{A}{\bar {x}}_{A}+n_{B}{\bar {x}}_{B}}{n_{AB}}}}$Reutilizando la WelfordVarianceclase anterior, tenemos:

def merge ( a : WelfordVariance , b : WelfordVariance ) - > WelfordVariance : ab = WelfordVariance ( ) ab.count = a.count + b.count delta = b.mean - a.mean ab.mean = ( a.count * a.mean + b.count * b.mean ) / ab.count ab.M2 = a.M2 + b.M2 + delta ** 2 * a.count * b.count / ab.count return ab​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​# ejemplo: ab = merge ( a , b ) print ( ab . get_sample_variance ())

Este algoritmo permite dividir un conjunto de datos en varias partes, ejecutarlas en paralelo y luego combinar los resultados. Esto posibilita la paralelización de cualquier forma, incluyendo AVX , con GPU y clústeres de computadoras . El algoritmo también puede adaptarse para estadísticas de orden superior, así como para covarianza. ^{[ 3 ] [ 12 ]}

Supongamos que todas las operaciones de punto flotante utilizan 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 insesgada de la varianza poblacional es 30. Tanto el algoritmo ingenuo como el algoritmo de dos pasadas calculan estos valores correctamente.

A continuación, consideremos 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 correctamente esta estimación de la varianza, pero el algoritmo ingenuo devuelve 29,333333333333332 en lugar de 30.

Si bien esta pérdida de precisión puede ser tolerable y considerarse un defecto menor del algoritmo ingenuo, aumentar aún más el desplazamiento hace que el error sea catastrófico. Consideremos 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 ^{[ 12 ]} extiende las fórmulas de Chan para calcular el tercer y cuarto momento central , necesarios, por ejemplo, para 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í el${\displaystyle M_{k}}$son nuevamente las sumas de potencias de diferencias con respecto a la media${\textstyle \sum (x-{\overline {x}})^{k}}$, donación

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

${\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${\displaystyle \delta /n}$Solo se necesita una operación de división y, por lo tanto, las estadísticas de orden superior se pueden calcular con un coste incremental mínimo.

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

def curtosis_en_línea ( datos ): n = media = M2 = M3 = M4 = 0para x en 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 puedes calcular la varianza usando M2 y la asimetría usando M3. # Precaución: si todas las entradas son iguales, M2 será 0, lo que resultará en una división por cero. curtosis = ( n * M4 ) / ( M2 ** 2 ) - 3 return curtosis

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

Choi y Sweetman ^{[ 15 ]} ofrecen dos métodos alternativos para calcular la asimetría y la curtosis, cada uno de los cuales puede ahorrar requisitos sustanciales de memoria de computadora y tiempo de CPU en ciertas aplicaciones. El primer enfoque consiste en calcular los momentos estadísticos separando los datos en intervalos 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 de orden superior. Una ventaja es que los cálculos de los momentos estadísticos se pueden realizar con una precisión arbitraria, de modo que los cálculos 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 la manera convencional: el rango de valores potenciales se divide en intervalos y el número de ocurrencias dentro de cada intervalo se cuenta y se grafica de manera que el área de cada rectángulo sea igual a la porción de los valores de la muestra dentro de ese intervalo:

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

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

${\displaystyle m_{n}^{(h)}=\sum _{k=1}^{K}x_{k}^{n}H(x_{k})\,\Delta x_{k}={\frac {1}{A}}\sum _{k=1}^{K}x_{k}^{n}h(x_{k})\,\Delta x_{k}}$

${\displaystyle \theta _{n}^{(h)}=\sum _{k=1}^{K}{\Big (}x_{k}-m_{1}^{(h)}{\Big )}^{n}\,H(x_{k})\,\Delta x_{k}={\frac {1}{A}}\sum _{k=1}^{K}{\Big (}x_{k}-m_{1}^{(h)}{\Big )}^{n}h(x_{k})\,\Delta x_{k}}$

donde el superíndice${\displaystyle ^{(h)}}$indica que los momentos se calculan a partir del histograma. Para un ancho de intervalo constante${\displaystyle \Delta x_{k}=\Delta x}$Estas dos expresiones se pueden simplificar usando${\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, propuesto por Choi y Sweetman ^{[ 15 ],} es una metodología analítica para combinar momentos estadísticos de segmentos individuales de una serie temporal, de manera que los momentos globales resultantes correspondan a la serie temporal completa. Esta metodología podría utilizarse para el cálculo paralelo de momentos estadísticos con su posterior combinación, o para la combinación de momentos estadísticos calculados en instantes de tiempo consecutivos.

Si${\displaystyle Q}$Se conocen conjuntos de momentos estadísticos: ${\displaystyle (\gamma _{0,q},\mu _{q},\sigma _{q}^{2},\alpha _{3,q},\alpha _{4,q})\quad }$para${\displaystyle q=1,2,\ldots ,Q}$, luego cada uno${\displaystyle \gamma _{n}}$puede expresarse en términos del equivalente${\displaystyle n}$momentos crudos:

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

dónde${\displaystyle \gamma _{0,q}}$Generalmente se considera que es la duración de la${\displaystyle q^{th}}$historial temporal, o el número de puntos si${\displaystyle \Delta t}$es constante.

El beneficio de expresar los momentos estadísticos en términos de${\displaystyle \gamma }$es que el${\displaystyle Q}$Los conjuntos se pueden combinar mediante la suma, y ​​no hay límite superior en el valor de${\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${\displaystyle _{c}}$representa el historial temporal concatenado o combinado${\displaystyle \gamma }$Estos valores combinados de${\displaystyle \gamma }$Luego se pueden transformar inversamente en momentos brutos que representan la historia temporal concatenada completa.

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

Relaciones conocidas entre los momentos crudos (${\displaystyle m_{n}}$) y los momentos centrales (${\displaystyle \theta _{n}=\operatorname {E} [(x-\mu )^{n}])}$Luego se utilizan 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 \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 .

def naive_covariance ( data1 , data2 ): n = len ( data1 ) sum1 = sum ( data1 ) sum2 = sum ( data2 ) sum12 = sum ([ i1 * i2 for i1 , i2 in zip ( data1 , data2 )])covarianza = ( suma12 - suma1 * suma2 / n ) / n devolver covarianza

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 ): # para mayor precisión, use sum(): Ex += ix - kx # Ex = sum(ix - kx for ix in data_x) # o sum(ix) - kx * n Ey += iy - ky # Ey = sum(iy - ky for iy in data_y) # o sum(iy) - ky * n Exy += ( ix - kx ) * ( iy - ky ) # Exy = sum((ix - kx) * (iy - ky) for ix, iy in zip(data_x, data_y)) return ( Exy - Ex * Ey / n ) / n

def two_pass_covariance ( data1 , data2 ): n = len ( data1 ) media1 = suma ( data1 ) / n media2 = suma ( data2 ) / ncovarianza = 0 para i1 , i2 en zip ( data1 , data2 ): # para mayor precisión, usar sum(): a = i1 - media1 # covarianza = sum((i1 - media1) * (i2 - media2) para i1, i2 en zip(data1, data2)) b = i2 - media2 covarianza += a * b / n return covarianza

def covarianza_en_línea ( datos1 , datos2 ): media_x = media_y = C = n = 0 para x , y en zip ( datos1 , datos2 ): n += 1 dx = x - media_x media_x += dx / n media_y += ( y - media_y ) / n C += dx * ( y - media_y )# Covarianza y corrección de Bessel para el retorno de la muestra ( C / n , C / ( n - 1 ))

def covarianza_ponderada_en_línea ( datos1 , datos2 , datos3 ): media_x = media_y = 0 suma_w = suma_w2 = 0 C = 0 para x , y , w en zip ( datos1 , datos2 , datos3 ): suma_w += w suma_w2 += w * w dx = x - media_x media_x += ( w / suma_w ) * dx media_y += ( w / suma_w ) * ( y - media_y ) C += w * dx * ( y - media_y )covarianza_población = C / wsum # Corrección de Bessel para la varianza de la muestra # Pesos de frecuencia covarianza_frecuencia_muestra = C / ( wsum - 1 ) # Pesos de fiabilidad covarianza_fiabilidad_muestra = C / ( wsum - wsum2 / wsum )

• Algoritmo de suma de Kahan
• Desviaciones cuadráticas respecto a la media
• Método Yamartino

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