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?
statistics
real-time
measurement
jonsca
fonte
fonte
Respostas:
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:
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:
e depois de cada etapa, o valor de
m
é a média e o desvio padrão pode ser calculado comosqrt(S/n)
ousqrt(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:
resultado:
Você notará que ainda há algum erro de arredondamento, mas não é ruim, enquanto
naive_stats
apenas vomita.edit: Acabei de notar o comentário de Belisarius citando a Wikipedia, que menciona o algoritmo de Knuth.
fonte
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
mean
emeansq
sejam as estimativas atuais da média e média do quadrado do sinal. A cada ciclo, atualize essas estimativas com a nova amostrax
:Aqui é uma constante que determina o comprimento efetivo da média corrente. Como escolher é descrito abaixo em "análise".a0<a<1 a
O que é expresso acima como um programa imperativo também pode ser representado como um diagrama de fluxo de sinal:
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+(1−a)yi−1 xi i yi z
Condensando os filtros IIR em seus próprios blocos, o diagrama agora fica assim:
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=esT T fs=1/T 1−(1−a)e−sT=0 s = 1Tregistro( 1 - a )
Escolha :uma
Referências
fonte
0 > a > 1
0 < a < 1
. Se o seu sistema possui um sistema de amostragemT
e você deseja uma constante de tempo médiotau
, escolhaa = 1 - exp (2*pi*T/tau)
.z=1
(DC) emH(z) = a/(1-(1-a)/z)
e você recebe 1.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:
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:i i
ou você pode usar:
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 :
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.
fonte
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
mean
es
é atualizado (variáveis de membro "globais"), igualm
es
acima na postagem de Jason. O valorcount
refere-se ao tamanho da janelan
.fonte
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.
fonte
y1 = filter(a,[1 (1-a)],x);
.