xticks barplot python perl statistics

python - barplot - pandas plot



¿Cómo calcular eficientemente una desviación estándar corriente? (14)

Tengo una variedad 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 eficientemente 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 recorriendo el conjunto y sumando 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, vuelvo a pasar, ahora que tengo la media calculada.

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

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

¿Qué tan grande es tu matriz? A menos que sean tropecientos de elementos largos, no te preocupes por recorrerlos dos veces. El código es simple y fácil de probar.

Mi preferencia sería usar la extensión numpy array maths para convertir su matriz de matrices en una matriz numpy 2D 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 eso no es una opción y necesitas una solución pura de Python, sigue 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 una sola vez, las sumas 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 listas anterior.

Aquí hay un "un-liner", distribuido en múltiples líneas, en un estilo de programació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)))

Aquí hay una traducción literal de Python pura 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();

Eche un vistazo a PDL (pronunciado "piddle!").

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

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;


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, tal vez sea PRMS (que muestran las estadísticas de Sinan :: ejemplo descriptivo) o RMS (que muestra el ejemplo NumPy de ars). Creo que uno de estos tres debe tener razón ;-)

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

El módulo Python Runstats es solo para este tipo de cosas. Instale runstats de 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 estadísticos 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 el Art of Computer Programming, Vol 2, p. 232, 3ª edición. El beneficio de esto es numéricamente estable y resultados precisos.

Descargo de responsabilidad: soy el autor del módulo Python Runstats.

La respuesta básica es acumular la suma de ambas x (llámala ''suma_x1'') yx2 (llámala ''suma_x2'') sobre la marcha. 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; obtienes la desviación estándar de la población usando ''n'' en lugar de ''n - 1'' como el 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 más información.

La respuesta es usar el algoritmo de Welford, que está muy claramente definido después de los "métodos ingenuos" en:

Es más numéricamente estable que los coleccionistas de dos pasos o de suma simple de cuadrados sugeridos en otras respuestas. La estabilidad solo importa cuando tienes muchos valores que están cerca uno del otro 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 cuadrada). Dividir por N-1 conduce a una estimación imparcial de la varianza de la muestra, mientras que dividir por N en promedio subestima la varianza (porque no tiene en cuenta la varianza entre la media de la muestra y la media real).

Escribí dos entradas de blog sobre el tema que incluyen más detalles, que incluyen cómo eliminar valores anteriores en línea:

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

Me gusta expresar 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 vea 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 de la muestra (1 / N), no la estimación no sesgada de la varianza de la población (que utiliza un factor de normalización 1 / (N-1)). A diferencia de las otras respuestas, la variable, var , que está rastreando la varianza de 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 el momento (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)

Quizás no sea lo que estabas preguntando, pero ... Si usas una matriz numpy, hará el trabajo por ti, eficientemente:

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 entrada de blog y comentarios sobre métodos de una sola pasada para medios y variaciones de computación:

Statistics::Descriptive es un módulo Perl muy decente para este tipo 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

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