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

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 recorriendo la matriz y sumando el valor en un índice dado de una lista. Al final, divido cada valor en mi "lista de promedios" por n(estoy trabajando con una población, no con una muestra de la población).

Para hacer la desviación estándar, vuelvo a recorrer, 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 utilizar 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 colectores de dos pasadas o la suma simple de cuadrados en línea 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 coma flotante.

Es posible que también 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). Dividir por N-1 conduce a una estimación insesgada 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 verdadera).

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

• Calcular la media y la varianza de la muestra en línea en una sola pasada
• Eliminación de valores en el algoritmo de Welford para la media y la varianza en línea

También puede echar un vistazo a mi implementación de Java; las pruebas javadoc, fuente y unitaria están todas en línea:

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

La respuesta básica es acumular la suma de x (llámelo 'sum_x1') y x ² (llámelo 'sum_x2') a medida que avanza. El valor de la desviación estándar es entonces:

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

dónde

mean = sum_x / n

Ésta 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 divisor.

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

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

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

import math

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()

print(f'Mean: {mean}, Variance: {variance}, Std. Dev.: {stdev}')

Quizás 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 un solo paso para calcular medias y variaciones:

• http://lingpipe-blog.com/2009/03/19/computing-sample-mean-variance-online-one-pass/

El módulo runstats de Python es solo para este tipo de cosas. Instale 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 estadísticos se basan en el método de Knuth y Welford para calcular la desviación estándar en una pasada, como se describe en Art of Computer Programming, Vol 2, p. 232, 3ª edición. El beneficio de esto son resultados numéricamente estables y precisos.

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

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

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

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

Aquí hay un ejemplo usando sus cifras ...

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, puede ser PRMS (que muestra el ejemplo de Estadística :: Descriptivo de Sinan) o RMS (que muestra el ejemplo de NumPy de ars). Supongo que uno de estos tres debe tener razón ;-)

Para obtener más información sobre la PDL, consulte:

• 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
• Wiki de PDL
• Entrada de Wikipedia para PDL
• Página del proyecto de Sourceforge para PDL

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

Mi preferencia sería usar la extensión matemática de matrices 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 pura de Python, siga leyendo ...

Si su 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 sumas corrientes 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.

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

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

Creo que este tema te ayudará. Desviación Estándar

Aquí hay un "resumen", distribuido en varias líneas, en 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)))

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 acumulada o en ejecución . Para eso, tendrá que convertir sus datos en un marco de datos de pandas (o una serie si es 1D), pero hay funciones para eso.

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 una sola pasada 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 se trata de calcular la varianza de la muestra (1 / N), no la estimación insesgada 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, varque sigue la varianza corriente 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)

Aquí hay un ejemplo práctico de cómo podría implementar una desviación estándar en ejecución con python y numpy:

a = np.arange(1, 10)
s = 0
s2 = 0
for i in range(0, len(a)):
    s += a[i]
    s2 += a[i] ** 2 
    n = (i + 1)
    m = s / n
    std = np.sqrt((s2 / n) - (m * m))
    print(std, np.std(a[:i + 1]))

Esto imprimirá la desviación estándar calculada y una desviación estándar de verificación calculada con numpy:

0.0 0.0
0.5 0.5
0.8164965809277263 0.816496580927726
1.118033988749895 1.118033988749895
1.4142135623730951 1.4142135623730951
1.707825127659933 1.707825127659933
2.0 2.0
2.29128784747792 2.29128784747792
2.5819888974716116 2.581988897471611

Solo estoy usando la fórmula descrita en este hilo:

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