Como calcular com eficiência um desvio padrão em execução?

Eu tenho uma série de listas de números, por exemplo:

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

O que eu gostaria de fazer é calcular com eficiência a média e o desvio padrão em cada índice de uma lista, em todos os elementos da matriz.

Para fazer a média, estive percorrendo a matriz e somando o valor em um determinado índice de uma lista. No final, divido cada valor em minha "lista de médias" por n(estou trabalhando com uma população, não uma amostra da população).

Para fazer o desvio padrão, faço um loop novamente, agora que calculei a média.

Eu gostaria de evitar percorrer a matriz duas vezes, uma para a média e outra para o SD (depois de ter uma média).

Existe um método eficiente para calcular os dois valores, passando pelo array apenas uma vez? Qualquer código em uma linguagem interpretada (por exemplo, Perl ou Python) ou pseudocódigo está bem.

A resposta é usar o algoritmo de Welford, que é muito claramente definido após os "métodos ingênuos" em:

• Wikipedia: Algoritmos para calcular a variância

É mais estável numericamente do que os coletores de soma de quadrados simples de duas passagens ou online sugeridos em outras respostas. A estabilidade só importa realmente quando você tem muitos valores próximos uns dos outros, pois eles levam ao que é conhecido como " cancelamento catastrófico " na literatura de ponto flutuante.

Você também pode querer melhorar a diferença entre dividir pelo número de amostras (N) e N-1 no cálculo de variância (desvio ao quadrado). A divisão por N-1 leva a uma estimativa imparcial da variância da amostra, enquanto a divisão por N na média subestima a variância (porque não leva em consideração a variância entre a média da amostra e a média verdadeira).

Escrevi duas entradas de blog sobre o tópico que fornecem mais detalhes, incluindo como excluir valores anteriores online:

• Calculando a média e a variação da amostra online em uma passagem
• Excluindo Valores no Algoritmo de Welford para Média e Variância Online

Você também pode dar uma olhada no meu implemento Java; o javadoc, o código-fonte e os testes de unidade estão todos online:

• Javadoc: stats.OnlineNormalEstimator
• Fonte: stats.OnlineNormalEstimator.java
• Fonte JUnit: test.unit.stats.OnlineNormalEstimatorTest.java
• Página inicial do LingPipe

A resposta básica é acumular a soma de x (chame de 'soma_x1') e x ² (chame de 'soma_x2') conforme você avança. O valor do desvio padrão é então:

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

Onde

mean = sum_x / n

Este é o desvio padrão da amostra; você obtém o desvio padrão da população usando 'n' em vez de 'n - 1' como divisor.

Você pode precisar se preocupar com a estabilidade numérica de tirar a diferença entre dois números grandes se estiver lidando com amostras grandes. Vá para as referências externas em outras respostas (Wikipedia, etc) para obter mais informações.

Aqui está uma tradução literal pura do Python da implementação do 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}')

Talvez não seja o que você estava pedindo, mas ... Se você usar uma matriz numpy, ela fará o trabalho para você, de forma 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 ]

A propósito, há uma discussão interessante nesta postagem do blog e comentários sobre métodos de uma passagem para meios de computação e variações:

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

O módulo runstats do Python serve exatamente para esse tipo de coisa. Instale runstats do PyPI:

pip install runstats

Os resumos de runstats podem produzir a média, a variação, o desvio padrão, a assimetria e a curtose em uma única passagem de dados. Podemos usar isso para criar sua versão "em execução".

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

Os resumos de estatísticas são baseados no método de Knuth e Welford para calcular o desvio padrão em uma passagem, conforme descrito em Art of Computer Programming, Vol 2, p. 232, 3ª edição. O benefício disso são resultados numericamente estáveis ​​e precisos.

Isenção de responsabilidade: eu sou o autor do módulo runstats do Python.

Statistics :: Descriptive é um módulo Perl muito decente para estes 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__

Resultado:

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

Dê uma olhada no PDL (pronuncia-se "piddle!").

Esta é a linguagem de dados Perl projetada para matemática de alta precisão e computação científica.

Aqui está um exemplo usando suas 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 produz:

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]

Dê uma olhada em PDL :: Primitive para obter mais informações sobre o statsover função . Isso parece sugerir que ADEV é o "desvio padrão".

No entanto, pode ser PRMS (que mostra o exemplo de Estatísticas :: Descritivo de Sinan) ou RMS (que mostra o exemplo de NumPy de ars). Acho que um desses três deve estar certo ;-)

Para obter mais informações sobre o PDL, dê uma olhada em:

• pdl.perl.org (página oficial do PDL).
• Guia de referência rápida de PDL em PerlMonks
• Artigo do Dr. Dobb sobre PDL
• PDL Wiki
• Entrada da Wikipedia para PDL
• Página do projeto do Sourceforge para PDL

Qual é o tamanho do seu array? A menos que tenha zilhões de elementos, não se preocupe em percorrê-lo duas vezes. O código é simples e facilmente testado.

Minha preferência seria usar a extensão matemática numpy array para converter seu array de arrays em um array 2D numpy e obter o desvio padrão diretamente:

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

Se isso não for uma opção e você precisar de uma solução Python pura, continue lendo ...

Se o seu array é

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

Então, o desvio padrão é:

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

Se você estiver determinado a percorrer sua matriz apenas uma vez, as somas corridas podem ser combinadas.

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

Isso não é tão elegante quanto a solução de compreensão de lista acima.

Você pode consultar o artigo da Wikipedia sobre Desvio Padrão , em particular a seção sobre Métodos de cálculo rápido.

Também descobri um artigo que usa Python, você deve ser capaz de usar o código dele sem muitas alterações: Mensagens subliminares - executando desvios padrão .

Acho que esse problema vai te ajudar. Desvio padrão

Aqui está um "one-liner", espalhado por várias linhas, em estilo de programação 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 descreve a seguinte resposta: O pandas / scipy / numpy fornece uma função de desvio padrão cumulativo? O módulo Python Pandas contém um método para calcular o desvio padrão em execução ou cumulativo . Para isso, você terá que converter seus dados em um dataframe do pandas (ou em uma série se for 1D), mas existem funções para isso.

Gosto de expressar a atualização desta forma:

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 uma função de passagem única se pareça com isto:

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)

observe que isso está calculando a variância da amostra (1 / N), não a estimativa não enviesada da variância da população (que usa um fator de normalização 1 / (N-1)). Ao contrário das outras respostas, a variável, varque acompanha a variância em execução, não cresce em proporção ao número de amostras. Em todos os momentos, é apenas a variância do conjunto de amostras visto até agora (não há "divisão final por n" para obter a variância).

Em uma aula seria assim:

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

Isso também funciona para amostras 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)