¿Cuál sería la forma ideal de encontrar la media y la desviación estándar de una señal para una aplicación en tiempo real? Me gustaría poder activar un controlador cuando una señal estuvo a más de 3 desviaciones estándar de la media durante un cierto período de tiempo.
Supongo que un DSP dedicado haría esto con bastante facilidad, pero ¿hay algún "atajo" que no requiera algo tan complicado?
statistics
real-time
measurement
jonsca
fuente
fuente
Respuestas:
Hay una falla en la respuesta de Jason R, que se discute en el vol. "Arte de la programación de computadoras" de Knuth. 2. El problema surge si tiene una desviación estándar que es una pequeña fracción de la media: el cálculo de E (x ^ 2) - (E (x) ^ 2) sufre de una severa sensibilidad a los errores de redondeo de coma flotante.
Incluso puedes probar esto tú mismo en un script de Python:
Obtengo -128.0 como respuesta, lo que claramente no es computacionalmente válido, ya que las matemáticas predicen que el resultado no debería ser negativo.
Knuth cita un enfoque (no recuerdo el nombre del inventor) para calcular la media de carrera y la desviación estándar que es algo así:
y luego, después de cada paso, el valor de
m
es la media, y la desviación estándar se puede calcular comosqrt(S/n)
osqrt(S/n-1)
según cuál sea su definición favorita de desviación estándar.La ecuación que escribo arriba es ligeramente diferente a la de Knuth, pero es computacionalmente equivalente.
Cuando tenga unos minutos más, codificaré la fórmula anterior en Python y mostraré que obtendrás una respuesta no negativa (que con suerte está cerca del valor correcto).
actualización: aquí está.
test1.py:
resultado:
Notarás que todavía hay algún error de redondeo, pero no está mal, mientras que
naive_stats
solo vomita.editar: Acabo de notar el comentario de Belisario citando Wikipedia que menciona el algoritmo de Knuth.
fuente
En el dominio de la frecuencia, un "promedio de ejecución ponderado exponencialmente" es simplemente un polo real. Es simple de implementar en el dominio del tiempo.
Implementación de dominio de tiempo
Sea
mean
ymeansq
sea las estimaciones actuales de la media y la media del cuadrado de la señal. En cada ciclo, actualice estas estimaciones con la nueva muestrax
:Lo que se expresa arriba como un programa imperativo también se puede representar como un diagrama de flujo de señal:
Análisis
Condensando los filtros IIR en sus propios bloques, el diagrama ahora se ve así:
Referencias
fuente
0 > a > 1
0 < a < 1
. Si su sistema tiene tmie de muestreoT
y desea una constante de tiempo promediotau
, elijaa = 1 - exp (2*pi*T/tau)
.z=1
(DC) enH(z) = a/(1-(1-a)/z)
y obtendrá 1.Un método que he usado antes en una aplicación de procesamiento integrada es mantener acumuladores de la suma y la suma de los cuadrados de la señal de interés:
o puedes usar:
dependiendo del método de estimación de desviación estándar que prefiera . Estas ecuaciones se basan en la definición de la varianza :
Los he usado con éxito en el pasado (aunque solo me preocupaba la estimación de la varianza, no la desviación estándar), aunque debe tener cuidado con los tipos numéricos que usa para mantener los acumuladores si va a sumar un largo periodo de tiempo; No quieres desbordamiento.
Editar: además del comentario anterior sobre el desbordamiento, debe tenerse en cuenta que este no es un algoritmo numéricamente robusto cuando se implementa en aritmética de punto flotante, lo que puede causar grandes errores en las estadísticas estimadas. Mire la respuesta de Jason S para un mejor enfoque en ese caso.
fuente
Similar a la respuesta preferida anterior (Jason S.), y también derivada de la fórmula tomada de Knut (Vol.2, p 232), también se puede derivar una fórmula para reemplazar un valor, es decir, eliminar y agregar un valor en un solo paso . Según mis pruebas, el reemplazo ofrece una mejor precisión que la versión de quitar / agregar de dos pasos.
El siguiente código está en Java
mean
ys
se actualiza (variables miembro "globales"), igual quem
y máss
arriba en la publicación de Jason. El valor secount
refiere al tamaño de la ventanan
.fuente
La respuesta de Jason y Nibot difiere en un aspecto importante: el método de Jason calcula el desarrollo estándar y la media de toda la señal (ya que y = 0), mientras que el de Nibot es un cálculo "en ejecución", es decir, pesa muestras más recientes más fuertes que las muestras de pasado lejano.
Dado que la aplicación parece requerir std dev y mean en función del tiempo, el método de Nibot es probablemente el más apropiado (para esta aplicación específica). Sin embargo, la parte realmente difícil será obtener la parte de ponderación de tiempo correcta. El ejemplo de Nibot usa un filtro simple de un polo.
La elección del filtro de paso bajo puede guiarse por lo que sabe sobre su señal y la resolución de tiempo que necesita para su estimación. Las frecuencias de corte más bajas y el orden más alto resultarán en una mejor precisión pero un tiempo de respuesta más lento.
Para complicar más las cosas, se aplica un filtro en el dominio lineal y otro en el dominio cuadrado. La cuadratura cambia significativamente el contenido espectral de la señal, por lo que es posible que desee utilizar un filtro diferente en el dominio cuadrado.
Aquí hay un ejemplo sobre cómo estimar la media, rms y std dev en función del tiempo.
fuente
y1 = filter(a,[1 (1-a)],x);
.