¿Cómo calcular de manera eficiente una desviación estándar de funcionamiento?

Tengo una serie de listas de números, por ejemplo:

[0] (0.01, 0.01, 0.02, 0.04, 0.03) [1] (0.00, 0.02, 0.02, 0.03, 0.02) [2] (0.01, 0.02, 0.02, 0.03, 0.02) ... [n] (0.01, 0.00, 0.01, 0.05, 0.03) 

Lo que me gustaría hacer es calcular de manera eficiente la media y la desviación estándar en cada índice de una lista, en todos los elementos de la matriz.

Para hacer la media, he estado haciendo un bucle a través de la matriz y sumndo el valor en un índice dado de una lista. Al final, divido cada valor en mi “lista de promedios” por n .

Para hacer la desviación estándar, repito, ahora que tengo la media calculada.

Me gustaría evitar pasar por la matriz dos veces, una para la media y luego una para la SD (después de tener una media).

¿Existe un método eficiente para calcular ambos valores, solo pasando por la matriz una vez? Cualquier código en un lenguaje interpretado (por ejemplo, Perl o Python) o pseudocódigo está bien.

La respuesta es usar el algoritmo de Welford, que se define muy claramente después de los “métodos ingenuos” en:

  • Wikipedia: Algoritmos para calcular la varianza.

Es más estable numéricamente que los recolectores de sum de cuadrados simples o en línea de dos pases sugeridos en otras respuestas. La estabilidad solo importa realmente cuando tienes muchos valores que están cerca unos de otros, ya que conducen a lo que se conoce como ” cancelación catastrófica ” en la literatura de punto flotante.

También es posible que desee repasar la diferencia entre dividir por el número de muestras (N) y N-1 en el cálculo de la varianza (desviación al cuadrado). La división por N-1 conduce a una estimación imparcial de la varianza de la muestra, mientras que al dividir por N en promedio se subestima la varianza (porque no tiene en cuenta la varianza entre la media de la muestra y la media verdadera).

Escribí dos entradas de blog sobre el tema que incluyen más detalles, incluida la forma de eliminar valores anteriores en línea:

  • Cálculo de la media y la varianza de la muestra en línea en un pase
  • Eliminar valores en el algoritmo de Welford para la media y la varianza en línea

También puedes echar un vistazo a mi implemento de Java; Las pruebas de javadoc, fuente y unidad están todas en línea:

  • Javadoc: stats.OnlineNormalEstimator
  • Fuente: stats.OnlineNormalEstimator.java
  • Fuente de JUnit: test.unit.stats.OnlineNormalEstimatorTest.java
  • Página de inicio de LingPipe

La respuesta básica es acumular la sum de ambos x (llámalo ‘sum_x1’) y x 2 (llámalo ‘sum_x2’) a medida que avanzas. El valor de la desviación estándar es entonces:

 stdev = sqrt((sum_x2 / n) - (mean * mean)) 

dónde

 mean = sum_x / n 

Esta es la desviación estándar de la muestra; se obtiene la desviación estándar de la población utilizando ‘n’ en lugar de ‘n – 1’ como divisor.

Es posible que deba preocuparse por la estabilidad numérica de tomar la diferencia entre dos números grandes si está tratando con muestras grandes. Vaya a las referencias externas en otras respuestas (Wikipedia, etc.) para obtener más información.

Tal vez no sea lo que estaba preguntando, pero … Si usa una matriz numpy, hará el trabajo por usted de manera eficiente:

 from numpy import array nums = array(((0.01, 0.01, 0.02, 0.04, 0.03), (0.00, 0.02, 0.02, 0.03, 0.02), (0.01, 0.02, 0.02, 0.03, 0.02), (0.01, 0.00, 0.01, 0.05, 0.03))) print nums.std(axis=1) # [ 0.0116619 0.00979796 0.00632456 0.01788854] print nums.mean(axis=1) # [ 0.022 0.018 0.02 0.02 ] 

Por cierto, hay una discusión interesante en esta publicación de blog y comentarios sobre métodos de una pasada para calcular medios y variaciones:

Aquí hay una traducción literal en Python de la implementación del algoritmo de Welford desde http://www.johndcook.com/standard_deviation.html :

https://github.com/liyanage/python-modules/blob/master/running_stats.py

 class RunningStats: def __init__(self): self.n = 0 self.old_m = 0 self.new_m = 0 self.old_s = 0 self.new_s = 0 def clear(self): self.n = 0 def push(self, x): self.n += 1 if self.n == 1: self.old_m = self.new_m = x self.old_s = 0 else: self.new_m = self.old_m + (x - self.old_m) / self.n self.new_s = self.old_s + (x - self.old_m) * (x - self.new_m) self.old_m = self.new_m self.old_s = self.new_s def mean(self): return self.new_m if self.n else 0.0 def variance(self): return self.new_s / (self.n - 1) if self.n > 1 else 0.0 def standard_deviation(self): return math.sqrt(self.variance()) 

Uso:

 rs = RunningStats() rs.push(17.0); rs.push(19.0); rs.push(24.0); mean = rs.mean(); variance = rs.variance(); stdev = rs.standard_deviation(); 

El módulo runstats de Python es solo para este tipo de cosas. Instala runstats desde PyPI:

 pip install runstats 

Los resúmenes de Runstats pueden producir la media, la varianza, la desviación estándar, la asimetría y la curtosis en una sola pasada de datos. Podemos usar esto para crear su versión “en ejecución”.

 from runstats import Statistics stats = [Statistics() for num in range(len(data[0]))] for row in data: for index, val in enumerate(row): stats[index].push(val) for index, stat in enumerate(stats): print 'Index', index, 'mean:', stat.mean() print 'Index', index, 'standard deviation:', stat.stddev() 

Los resúmenes de estadísticas se basan en el método de Knuth y Welford para calcular la desviación estándar en una sola pasada, tal como se describe en Art of Computer Programming, Vol. 2, pág. 232, 3ª edición. El beneficio de esto es resultados numéricamente estables y precisos.

Descargo de responsabilidad: Soy el autor del módulo de estadísticas de Python.

Echa un vistazo a PDL (pronunciado “piddle!”).

Este es el lenguaje de datos Perl, diseñado para matemáticas de alta precisión y computación científica.

Aquí hay un ejemplo usando tus figuras …

 use strict; use warnings; use PDL; my $figs = pdl [ [0.01, 0.01, 0.02, 0.04, 0.03], [0.00, 0.02, 0.02, 0.03, 0.02], [0.01, 0.02, 0.02, 0.03, 0.02], [0.01, 0.00, 0.01, 0.05, 0.03], ]; my ( $mean, $prms, $median, $min, $max, $adev, $rms ) = statsover( $figs ); say "Mean scores: ", $mean; say "Std dev? (adev): ", $adev; say "Std dev? (prms): ", $prms; say "Std dev? (rms): ", $rms; 

Lo que produce:

 Mean scores: [0.022 0.018 0.02 0.02] Std dev? (adev): [0.0104 0.0072 0.004 0.016] Std dev? (prms): [0.013038405 0.010954451 0.0070710678 0.02] Std dev? (rms): [0.011661904 0.009797959 0.0063245553 0.017888544] 

Eche un vistazo a PDL :: Primitive para obtener más información sobre la función statsover . Esto parece sugerir que ADEV es la “desviación estándar”.

Sin embargo, puede que sea PRMS (que muestra el ejemplo descriptivo de Sinan Statistics ::) o RMS (que muestra el ejemplo NumPy de ars). Supongo que uno de estos tres debe tener razón 😉

Para más información de PDL, eche un vistazo a:

  • pdl.perl.org (página oficial de PDL).
  • Guía de referencia rápida de PDL en PerlMonks
  • Artículo del Dr. Dobb sobre PDL
  • PDL Wiki
  • Entrada de Wikipedia para PDL
  • Página de proyecto de Sourceforge para PDL

Statistics :: Descriptivo es un módulo de Perl muy decente para estos tipos de cálculos:

 #!/usr/bin/perl use strict; use warnings; use Statistics::Descriptive qw( :all ); my $data = [ [ 0.01, 0.01, 0.02, 0.04, 0.03 ], [ 0.00, 0.02, 0.02, 0.03, 0.02 ], [ 0.01, 0.02, 0.02, 0.03, 0.02 ], [ 0.01, 0.00, 0.01, 0.05, 0.03 ], ]; my $stat = Statistics::Descriptive::Full->new; # You also have the option of using sparse data structures for my $ref ( @$data ) { $stat->add_data( @$ref ); printf "Running mean: %f\n", $stat->mean; printf "Running stdev: %f\n", $stat->standard_deviation; } __END__ 

Salida:

 C:\Temp> g Running mean: 0.022000 Running stdev: 0.013038 Running mean: 0.020000 Running stdev: 0.011547 Running mean: 0.020000 Running stdev: 0.010000 Running mean: 0.020000 Running stdev: 0.012566 

¿Qué tan grande es tu matriz? A menos que se trate de miles de millones de elementos, no se preocupe por recorrerlos dos veces. El código es simple y fácil de probar.

Mi preferencia sería usar la extensión matemática de matriz numpy para convertir su matriz de matrices en una matriz 2D numpy y obtener la desviación estándar directamente:

 >>> x = [ [ 1, 2, 4, 3, 4, 5 ], [ 3, 4, 5, 6, 7, 8 ] ] * 10 >>> import numpy >>> a = numpy.array(x) >>> a.std(axis=0) array([ 1. , 1. , 0.5, 1.5, 1.5, 1.5]) >>> a.mean(axis=0) array([ 2. , 3. , 4.5, 4.5, 5.5, 6.5]) 

Si esa no es una opción y necesita una solución de Python pura, siga leyendo …

Si tu matriz es

 x = [ [ 1, 2, 4, 3, 4, 5 ], [ 3, 4, 5, 6, 7, 8 ], .... ] 

Entonces la desviación estándar es:

 d = len(x[0]) n = len(x) sum_x = [ sum(v[i] for v in x) for i in range(d) ] sum_x2 = [ sum(v[i]**2 for v in x) for i in range(d) ] std_dev = [ sqrt((sx2 - sx**2)/N) for sx, sx2 in zip(sum_x, sum_x2) ] 

Si está decidido a recorrer su matriz solo una vez, las sums streams se pueden combinar.

 sum_x = [ 0 ] * d sum_x2 = [ 0 ] * d for v in x: for i, t in enumerate(v): sum_x[i] += t sum_x2[i] += t**2 

Esto no es tan elegante como la solución de comprensión de lista anterior.

Creo que este problema te ayudará. Desviación estándar

Puede consultar el artículo de Wikipedia sobre Desviación estándar , en particular la sección sobre Métodos de cálculo rápidos.

También hay un artículo que encontré que usa Python, debería poder usar el código sin muchos cambios: Mensajes subliminales : desviaciones estándar en ejecución .

 n=int(raw_input("Enter no. of terms:")) L=[] for i in range (1,n+1): x=float(raw_input("Enter term:")) L.append(x) sum=0 for i in range(n): sum=sum+L[i] avg=sum/n sumdev=0 for j in range(n): sumdev=sumdev+(L[j]-avg)**2 dev=(sumdev/n)**0.5 print "Standard deviation is", dev 

Como describe la siguiente respuesta: ¿Pandas / scipy / numpy proporciona una función de desviación estándar acumulativa? El módulo Python Pandas contiene un método para calcular la desviación estándar en ejecución o acumulada . Para eso, tendrás que convertir tus datos en un dataframe de pandas (o una serie si es 1D), pero hay funciones para eso.

Aquí hay un “one-liner”, distribuido en múltiples líneas, en un estilo de progtwigción funcional:

 def variance(data, opt=0): return (lambda (m2, i, _): m2 / (opt + i - 1))( reduce( lambda (m2, i, avg), x: ( m2 + (x - avg) ** 2 * i / (i + 1), i + 1, avg + (x - avg) / (i + 1) ), data, (0, 0, 0))) 

Me gusta express la actualización de esta manera:

 def running_update(x, N, mu, var): ''' @arg x: the current data sample @arg N : the number of previous samples @arg mu: the mean of the previous samples @arg var : the variance over the previous samples @retval (N+1, mu', var') -- updated mean, variance and count ''' N = N + 1 rho = 1.0/N d = x - mu mu += rho*d var += rho*((1-rho)*d**2 - var) return (N, mu, var) 

para que una función de un paso se vería así:

 def one_pass(data): N = 0 mu = 0.0 var = 0.0 for x in data: N = N + 1 rho = 1.0/N d = x - mu mu += rho*d var += rho*((1-rho)*d**2 - var) # could yield here if you want partial results return (N, mu, var) 

tenga en cuenta que esto es calcular la varianza muestral (1 / N), no la estimación imparcial de la varianza poblacional (que utiliza un factor de normalización 1 / (N-1)). A diferencia de las otras respuestas, la variable, var , que está siguiendo la varianza en ejecución no crece en proporción al número de muestras. En todo momento es solo la varianza del conjunto de muestras vistas hasta ahora (no hay una “división por n” final para obtener la varianza).

En una clase se vería así:

 class RunningMeanVar(object): def __init__(self): self.N = 0 self.mu = 0.0 self.var = 0.0 def push(self, x): self.N = self.N + 1 rho = 1.0/N d = x-self.mu self.mu += rho*d self.var += + rho*((1-rho)*d**2-self.var) # reset, accessors etc. can be setup as you see fit 

Esto también funciona para muestras ponderadas:

 def running_update(w, x, N, mu, var): ''' @arg w: the weight of the current sample @arg x: the current data sample @arg mu: the mean of the previous N sample @arg var : the variance over the previous N samples @arg N : the number of previous samples @retval (N+w, mu', var') -- updated mean, variance and count ''' N = N + w rho = w/N d = x - mu mu += rho*d var += rho*((1-rho)*d**2 - var) return (N, mu, var)