¿Extraer puntos de datos de la media móvil?

15

¿Es posible extraer puntos de datos de datos de promedio móvil?

En otras palabras, si un conjunto de datos solo tiene promedios móviles simples de los 30 puntos anteriores, ¿es posible extraer los puntos de datos originales?

¿Si es así, cómo?


fuente
1
La respuesta es un sí calificado, pero el procedimiento exacto depende de cómo se trate el segmento inicial de datos. Si simplemente se descarta, efectivamente ha perdido 15 datos, dejándolo con un sistema indeterminado de ecuaciones lineales. El resultado es que existen muchas respuestas válidas en general, pero aún puede avanzar si (a) se usan ventanas más cortas (o algún procedimiento de este tipo) para los 15 promedios móviles iniciales o (b) puede especificar restricciones adicionales en la solución (aproximadamente 15 dimensiones de restricciones ...). ¿En qué situación estás?
whuber
@whuber Muchas gracias por mirar! Tengo 2,000 puntos. El primer punto MA probablemente sea un promedio de los primeros 30 puntos originales. La precisión es la segunda a un resultado generalmente correcto, más específicamente buenas conjeturas en los puntos más "recientes". ¿Me puede recomendar un método relativamente simple? ¡Gracias por adelantado!
1
(si toma más de cinco minutos para escribir un comentario ...). Lo que quería escribir es que puedes pensar en el promedio como una multiplicación matricial. Las filas en el medio tendrán 1/30 * [1 1 1 ...] antes de la diagonal. La pregunta es, ¿cómo manejas los puntos en los bordes de tu vector para hacer que la matriz sea invertible? Puede hacer esto asumiendo que son el resultado de promediar menos elementos o si piensa en otras restricciones. Tenga en cuenta que si bien una inversión matricial es una forma fácil de entenderla, no es la más eficiente. Probablemente quieras usar un FFT para hacer eso.
fabee

Respuestas:

4

+1 a la respuesta de fabee, que está completa. Solo una nota para traducirlo a R, en función de los paquetes que he encontrado para realizar las operaciones en cuestión. En mi caso, tenía datos que son pronósticos de temperatura NOAA por tres meses: enero-febrero-marzo, febrero-marzo-abril, marzo-abril-mayo, etc., y quería dividirlos en (aproximadamente) valores mensuales, suponiendo que la temperatura de cada período de tres meses es esencialmente un promedio.

library (Matrix)
library (matrixcalc)

# Feb-Mar-Apr through Nov-Dec-Jan temperature forecasts:

qtemps <- c(46.0, 56.4, 65.8, 73.4, 77.4, 76.2, 69.5, 60.1, 49.5, 41.2)

# Thus I need a 10x12 matrix, which is a band matrix but with the first
# and last rows removed so that each row contains 3 1's, for three months.
# Yeah, the as.matrix and all is a bit obfuscated, but the results of
# band are not what svd.inverse wants.

a <- as.matrix (band (matrix (1, nrow=12, ncol=12), -1, 1)[-c(1, 12),])
ai <- svd.inverse (a)

mtemps <- t(qtemps) %*% t(ai) * 3

Lo cual funciona muy bien para mí. Gracias @fabee.

EDITAR: OK, traduciendo mi R a Python, obtengo:

from numpy import *
from numpy.linalg import *

qtemps = transpose ([[46.0, 56.4, 65.8, 73.4, 77.4, 76.2, 69.5, 60.1, 49.5, 41.2]])

a = tril (ones ((12, 12)), 2) - tril (ones ((12, 12)), -1)
a = a[0:10,:]

ai = pinv (a)

mtemps = dot (ai, qtemps) * 3

(Lo que llevó mucho más tiempo depurar que la versión R. Primero porque no estoy tan familiarizado con Python como con R, pero también porque R es mucho más utilizable de forma interactiva).

Wayne
fuente
@Gracchus: Lo siento, no soy un chico de C ++, pero puedes encontrar lo que necesitas en la biblioteca de álgebra lineal Armadillo C ++ ( arma.sourceforge.net ), que también está disponible en R a través del paquete RcppArmadillo.
Wayne
OK, mira si te funciona. Si es así, puedes elegir mi respuesta ;-)
Wayne
Las mejores prácticas para su información en Python son hacer importaciones absolutas: python.org/dev/peps/pep-0008/#imports, lo que hace que sea mucho más fácil leer el código de otras personas, porque realmente sabe de dónde provienen las funciones en lugar de tener que mira cada uno que no conoces. Ojalá fuera estándar en R hacer lo mismo. Tener que buscar todas las pequeñas funciones en el código de otra persona realmente hace que mis engranajes ...
palabras para
Además, los cuadernos Jupyter para la interactividad de Python o IPython.
wordsforthewise
17

xn=2000=30y=Axx

A=130(1...10...001...10...0...1...100...01...1)

3030y19702000

x1,...,x2000y1y2

x1,...,xnxyx

A3030AA

AAz=AyxyAz

2000x

Reconstrucción de la señal original a partir de la media móvil utilizando el pseudoinverso

Muchos programas numéricos ofrecen pseudo-inversas (por ejemplo, Matlab, numpy en python, etc.).

Aquí estaría el código de Python para generar las señales de mi ejemplo:

from numpy import *
from numpy.linalg import *
from matplotlib.pyplot import *
# get A and its inverse     
A = (tril(ones((2000,2000)),-1) - tril(ones((2000,2000)),-31))/30.
A = A[30:,:]
pA = pinv(A) #pseudo inverse

# get x
x = random.randn(2000) + 5
y = dot(A,x)

# reconstruct
x2 = dot(pA,y)

plot(x,label='original x')
plot(y,label='averaged x')
plot(x2,label='reconstructed x')
legend()
show()

Espero que ayude.

fabee
fuente
Esta es una gran respuesta, pero creo que está equivocado cuando dijo que "minimiza la distancia cuadrática entre y y Az". De hecho, y y Az son lo mismo. Lo que se minimiza es la norma de z que funciona bien para las señales del mundo real que he probado, pero no es tan buena si su señal original tiene muchos valores atípicos.
gdelfino
No estoy seguro de seguirlo. y y Ax son lo mismo, pero no y y Az Es cierto que también minimiza la norma de z. Tampoco veo por qué no funciona para mis ejemplos. La línea azul y la roja coinciden bastante bien. ¿Me estoy perdiendo algo en tu comentario?
fabee
y es el promedio móvil calculado a partir de la señal original x multiplicando por A. Este procedimiento nos da una señal z que tiene el mismo promedio móvil y. Por lo tanto, y = Az Entonces, solo la norma de z se minimiza. Si la señal original tiene un gran valor normal, entonces el procedimiento no dará buenos resultados. A continuación se muestra una señal de ejemplo con un valor de norma grande:
gdelfino
{42.8, -33.7, 13.2, -45.6, 10.2, 35.8, -41.4, 20.253, 43.3429, -33.2735, 13.6135, -45.1067, 10.6346, 36.1352, -40.9703, 20.6616, 43.6796, -32.8966, 14.0406, -44.7001, 10.9988 , 36.4675, -40.7277, 20.8823, 43.7878, -32.7415, 13.9951, -44.7947, 11.044, 36.3873, -40.7117, 20.7505, 43.8204, -32.9399, 13.9129, -44.9549, 10.8703, 36.1559, -40.889, 43.478.42.41.459, -40.889. , 13.5468, -45.2374, 10.3787, 35.8235, -41.5161, 19.9717, 43.0658, -33.7125, 13.0321}
gdelfino
Utilice un tamaño de ventana de 8 para la señal anterior. De esta forma, la señal filtrada tiene una forma muy diferente de la señal original.
gdelfino