Determinando a média e o desvio padrão em tempo real

31

Qual seria a maneira ideal de encontrar a média e o desvio padrão de um sinal para uma aplicação em tempo real. Eu gostaria de poder acionar um controlador quando um sinal estivesse a mais de 3 desvios-padrão da média por um certo período de tempo.

Estou assumindo que um DSP dedicado faria isso facilmente, mas existe algum "atalho" que pode não exigir algo tão complicado?

jonsca
fonte
Você sabe alguma coisa sobre o sinal? É estacionário?
@ Tim Vamos dizer que é estacionário. Para minha própria curiosidade, quais seriam as ramificações de um sinal não estacionário?
amigos estão dizendo sobre jonsca
3
Se estiver parado, você pode simplesmente calcular uma média e um desvio padrão em execução. As coisas seriam mais complicadas se a média e o desvio padrão variassem com o tempo.
5
Muito relacionado: en.wikipedia.org/wiki/…
Dr. belisarius

Respostas:

36

Há uma falha na resposta de Jason R, que é discutida no "Art of Computer Programming", de Knuth, vol. 2. O problema surge se você tiver um desvio padrão que é uma pequena fração da média: o cálculo de E (x ^ 2) - (E (x) ^ 2) sofre de severa sensibilidade a erros de arredondamento de ponto flutuante.

Você pode até tentar isso sozinho em um script Python:

ofs = 1e9
A = [ofs+x for x in [1,-1,2,3,0,4.02,5]] 
A2 = [x*x for x in A]
(sum(A2)/len(A))-(sum(A)/len(A))**2

Recebo -128.0 como resposta, que claramente não é computacionalmente válida, pois a matemática prevê que o resultado deve ser não-negativo.

Knuth cita uma abordagem (não me lembro do nome do inventor) para calcular a média e o desvio padrão em execução, que são mais ou menos assim:

 initialize:
    m = 0;
    S = 0;
    n = 0;

 for each incoming sample x:
    prev_mean = m;
    n = n + 1;
    m = m + (x-m)/n;
    S = S + (x-m)*(x-prev_mean);

e depois de cada etapa, o valor de mé a média e o desvio padrão pode ser calculado como sqrt(S/n)ou sqrt(S/n-1)dependendo da sua definição favorita de desvio padrão.

A equação que escrevo acima é um pouco diferente da de Knuth, mas é computacionalmente equivalente.

Quando tiver mais alguns minutos, codificarei a fórmula acima em Python e mostrarei que você obterá uma resposta não negativa (que, esperamos, esteja próxima do valor correto).


atualização: aqui está.

test1.py:

import math

def stats(x):
  n = 0
  S = 0.0
  m = 0.0
  for x_i in x:
    n = n + 1
    m_prev = m
    m = m + (x_i - m) / n
    S = S + (x_i - m) * (x_i - m_prev)
  return {'mean': m, 'variance': S/n}

def naive_stats(x):
  S1 = sum(x)
  n = len(x)
  S2 = sum([x_i**2 for x_i in x])
  return {'mean': S1/n, 'variance': (S2/n - (S1/n)**2) }

x1 = [1,-1,2,3,0,4.02,5] 
x2 = [x+1e9 for x in x1]

print "naive_stats:"
print naive_stats(x1)
print naive_stats(x2)

print "stats:"
print stats(x1)
print stats(x2)

resultado:

naive_stats:
{'variance': 4.0114775510204073, 'mean': 2.0028571428571427}
{'variance': -128.0, 'mean': 1000000002.0028572}
stats:
{'variance': 4.0114775510204073, 'mean': 2.0028571428571431}
{'variance': 4.0114775868357446, 'mean': 1000000002.0028571}

Você notará que ainda há algum erro de arredondamento, mas não é ruim, enquanto naive_statsapenas vomita.


edit: Acabei de notar o comentário de Belisarius citando a Wikipedia, que menciona o algoritmo de Knuth.

Jason S
fonte
1
+1 para a resposta detalhada com código de exemplo. Essa abordagem é superior à indicada na minha resposta quando uma implementação de ponto flutuante é necessária.
Jason R
1
Pode-se também verificar isso para a implementação de um C ++: johndcook.com/standard_deviation.html
Rui Marques
1
sim, é isso. Ele usa as equações exatas que Knuth usa. Você pode otimizar um pouco e evitar a necessidade de verificar a iteração inicial x as iterações subsequentes, se você usar o meu método.
Jason S
"Knuth cita uma abordagem (não me lembro do nome do inventor) para calcular a média corrente" - a propósito, é o método de Welford .
Jason S
Eu postei uma pergunta relacionada a isso, se alguém puder ajudar: dsp.stackexchange.com/questions/31812/…
Jonathan
13

Qual seria a maneira ideal de encontrar a média e o desvio padrão de um sinal para uma aplicação em tempo real. Eu gostaria de poder acionar um controlador quando um sinal estivesse a mais de 3 desvios-padrão da média por um determinado período de tempo.

A abordagem correta em situações como essa geralmente é calcular uma média de execução ponderada exponencialmente e um desvio padrão. Na média ponderada exponencialmente, as estimativas da média e da variância são enviesadas para a amostra mais recente, fornecendo estimativas da média e da variância nos últimos segundosτ , que provavelmente é o que você deseja, em vez da média aritmética usual em todos os amostras já vistas.

No domínio da frequência, uma "média de corrida ponderada exponencialmente" é simplesmente um polo real. É simples de implementar no domínio do tempo.

Implementação no domínio do tempo

Seja meane meansqsejam as estimativas atuais da média e média do quadrado do sinal. A cada ciclo, atualize essas estimativas com a nova amostra x:

% update the estimate of the mean and the mean square:
mean = (1-a)*mean + a*x
meansq = (1-a)*meansq + a*(x^2)

% calculate the estimate of the variance:
var = meansq - mean^2;

% and, if you want standard deviation:
std = sqrt(var);

Aqui é uma constante que determina o comprimento efetivo da média corrente. Como escolher é descrito abaixo em "análise".a0<a<1a

O que é expresso acima como um programa imperativo também pode ser representado como um diagrama de fluxo de sinal:

insira a descrição da imagem aqui

Análise

O algoritmo acima calcula onde é a entrada na amostra e é a saída (ou seja, estimativa da média). Este é um filtro IIR simples, de um pólo. Tomando a transformação , encontramos a função de transferência . x i i y i z H ( z ) = ayi=axi+(1a)yi1xiiyiz

H(z)=a1(1a)z1

Condensando os filtros IIR em seus próprios blocos, o diagrama agora fica assim:

insira a descrição da imagem aqui

Para ir para o domínio contínuo, fazemos a substituição que é o tempo da amostra e é a taxa da amostra. Resolvendo , descobrimos que o sistema contínuo possui um polo em . z=esTTfs=1/T1(1a)esT=0s=1Tlog(1a)

Escolha :a

a=1exp{2πTτ}

Referências

nibot
fonte
1
Você poderia explicar como determina a duração da média atual? E que valor de deve ser usado? A especificação é impossível de atender. aa0 > a > 1
usar o seguinte código
Isso é semelhante à abordagem de Jason R. Este método será menos preciso, mas um pouco mais rápido e com menos memória. Essa abordagem acaba usando uma janela exponencial.
Schnarf
Woops! Claro que eu quis dizer 0 < a < 1. Se o seu sistema possui um sistema de amostragem Te você deseja uma constante de tempo médio tau, escolha a = 1 - exp (2*pi*T/tau).
Nibot
Eu acho que pode haver um erro aqui. Os filtros unipolares não têm ganho de 0 dB em CC e, como você aplica um filtro no domínio linear e outro no domínio ao quadrado, o erro de ganho é diferente para E <x> e E <x ^ 2>. Vou elaborar mais na minha resposta
Hilmar
Possui ganho de 0 dB no DC. Substitua z=1(DC) em H(z) = a/(1-(1-a)/z)e você recebe 1.
nibot
5

Um método que eu usei antes em um aplicativo de processamento incorporado é manter acumuladores da soma e soma dos quadrados do sinal de interesse:

Ax,i=k=0ix[k]=Ax,i1+x[i],Ax,1=0

Ax2,i=k=0ix2[k]=Ax2,i1+x2[i],Ax2,1=0

Além disso, acompanhe o instante de tempo atual nas equações acima (ou seja, anote o número de amostras que você adicionou nos acumuladores). Então, a média da amostra e o desvio padrão no momento são:ii

μ~=Axii+1

σ~=Axi2i+1μ~2

ou você pode usar:

σ~=Axi2iμ~2

dependendo de qual método de estimativa de desvio padrão você preferir . Estas equações são baseadas na definição da variância :

σ2=E(X2)(E(X))2

Eu as usei com sucesso no passado (embora eu estivesse preocupado apenas com a estimativa de variância, não com o desvio padrão), embora você tenha que ter cuidado com os tipos numéricos que usa para manter os acumuladores, se quiser fazer uma soma um longo período de tempo; você não quer transbordar.

Edit: Além do comentário acima sobre estouro, note-se que este não é um algoritmo numericamente robusto quando implementado na aritmética de ponto flutuante, causando potencialmente grandes erros nas estatísticas estimadas. Veja a resposta de Jason S para uma abordagem melhor nesse caso.

Jason R
fonte
1
Talvez reescreva-o como deixe claro que você precisa adicionar apenas dois números em cada etapa, para que ninguém a implemente como soma todos os elementos de em cada etapa. Ax,i=x[i]+Ax,i1, Ax,0=x[0]ix
Lorem Ipsum
Sim, está melhor. Tentei reescrever para tornar a implementação recursiva mais clara.
Jason R
2
-1 quando tenho representante suficiente para fazê-lo: isso tem problemas numéricos. Veja Knuth vol. 2
Jason S
Parece haver alguns erros de digitação aqui. Por que a média está sendo subtraída sob o sinal de raiz quadrada de ? deve ser para corresponder à equação exibida , não? Além disso, apesar de não votar nesta resposta, concordo com Jason S de que pode haver problemas numéricos nessa abordagem. σμ2σ2=E(X2)(E(X))2
usar o seguinte código
2
@ Jason: Eu discordo que a técnica é inerentemente falha, embora eu concorde com o seu ponto de vista que não é um método numericamente robusto quando implementado em ponto flutuante. Eu deveria ter deixado mais claro que usei isso com sucesso em um aplicativo que usava aritmética inteira . A aritmética inteira (ou implementações de ponto fixo de números fracionários) não sofre com o problema que você apontou que causa perda de precisão. Nesse contexto, é um método adequado que requer menos operações por amostra.
Jason R
3

Semelhante à resposta preferida acima (Jason S.) e também derivada da fórmula de Knut (Vol.2, p 232), também é possível derivar uma fórmula para substituir um valor, ou seja, remover e adicionar um valor em uma etapa . De acordo com meus testes, o replace fornece uma precisão melhor do que a versão remover / adicionar em duas etapas.

O código abaixo está em Java meane sé atualizado (variáveis ​​de membro "globais"), igual me sacima na postagem de Jason. O valor countrefere-se ao tamanho da janela n.

/**
 * Replaces the value {@code x} currently present in this sample with the
 * new value {@code y}. In a sliding window, {@code x} is the value that
 * drops out and {@code y} is the new value entering the window. The sample
 * count remains constant with this operation.
 * 
 * @param x
 *            the value to remove
 * @param y
 *            the value to add
 */
public void replace(double x, double y) {
    final double deltaYX = y - x;
    final double deltaX = x - mean;
    final double deltaY = y - mean;
    mean = mean + deltaYX / count;
    final double deltaYp = y - mean;
    final double countMinus1 = count - 1;
    s = s - count / countMinus1 * (deltaX * deltaX - deltaY * deltaYp) - deltaYX * deltaYp / countMinus1;
}
marco
fonte
3

A resposta de Jason e Nibot difere em um aspecto importante: o método de Jason calcula o desvio padrão e a média para todo o sinal (desde y = 0), enquanto o de Nibot é um cálculo "contínuo", ou seja, pesa amostras mais recentes mais fortes que as amostras do passado distante.

Como o aplicativo parece exigir std dev e mean em função do tempo, o método de Nibot é provavelmente o mais apropriado (para esse aplicativo específico). No entanto, a parte mais complicada será acertar a parte de ponderação do tempo. O exemplo de Nibot usa um filtro simples de um pólo.

A maneira correta de descrever isso é obter uma estimativa de filtrando e uma estimativa para filtrando . Os filtros de estimativa são geralmente filtros de passa-baixo. Esses filtros devem ser dimensionados para obter ganho de 0dB a 0 Hz. Caso contrário, há um erro de ganho constante.x [ n ] E [ x 2 ] x [ n ] 2E[x]x[n]E[x2]x[n]2

A escolha do filtro passa-baixo pode ser orientada pelo que você sabe sobre o seu sinal e pela resolução de tempo necessária para sua estimativa. Freqüências de corte mais baixas e ordem mais alta resultarão em melhor precisão, mas menor tempo de resposta.

Para complicar ainda mais, um filtro é aplicado no domínio linear e outro no domínio ao quadrado. O quadrado altera significativamente o conteúdo espectral do sinal, portanto, você pode usar um filtro diferente no domínio do quadrado.

Aqui está um exemplo de como estimar média, rms e std dev em função do tempo.

%% example
fs = 44100; n = fs; % 44.1 kHz sample rate, 1 second
% signal: white noise plus a low frequency drift at 5 Hz)
x = randn(n,1) + sin(2*pi*(0:n-1)'*5/fs);
% mean estimation filter: since we are looking for effects in the 5 Hz range we use maybe a
% 25 Hz filter, 2nd order so it's not too sluggish
[b,a] = butter(2,25*2/fs);
xmeanEst = filter(b,a,x);
% now we estimate x^2, since most frequency double we use twice the bandwidth
[b,a] = butter(2,50*2/fs);
x2Est = filter(b,a,x.^2);
% std deviation estimate
xstd = sqrt(x2Est)-xmeanEst;
% and plot it
h = plot([x, xmeanEst sqrt(x2Est) xstd]);
grid on;
legend('x','E<x>','sqrt(E<x^2>)','Std dev');
set(h(2:4),'Linewidth',2);
Hilmar
fonte
1
O filtro na minha resposta corresponde y1 = filter(a,[1 (1-a)],x);.
Nibot 01/03
1
Bom argumento sobre a distinção entre as estatísticas em execução e as estatísticas da amostra geral. Minha implementação pode ser modificada para calcular estatísticas em execução, acumulando-se sobre uma janela em movimento, o que também pode ser feito com eficiência (a cada etapa, subtraia a amostra de tempo que acabou de sair da janela de cada acumulador).
Jason R
nibot, desculpe você está certo e eu estava errado. Corrigirei isso imediatamente
Hilmar
1
1 para sugerindo filtragem diferente para X e X ^ 2
nibot