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

Question

¿Es posible extraer puntos de datos de la media móvil?

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

Si es así, ¿cómo?

Answer 1

Trato de poner lo que dijo Whuber en una respuesta. Digamos que tienes un gran vector $\mathbf x$ avec $n=2000$ entradas. Si se calcula una media móvil con una ventana de longitud $\ell=30$ se puede escribir como una multiplicación matricial vectorial $\mathbf y = A\mathbf x$ del vector $\mathbf x$ con la matriz

$$A=\frac{1}{30}\left(\begin{array}{cccccc} 1 & ... & 1 & 0 & ... & 0\\ 0 & 1 & ... & 1 & 0 & ...\\ \vdots & & \ddots & & & \vdots\\ 0 & ... & 1 & ... & 1 & 0\\ 0 & ... & 0 & 1 & ... & 1 \end{array}\right)$$

que tiene $30$ que se van desplazando a medida que se avanza en las filas hasta que el $30$ los que llegan al final de la matriz. Aquí el vector promediado $\mathbf y$ tiene las dimensiones de 1970. La matriz tiene $1970$ filas y $2000$ columnas. Por lo tanto, no es invertible.

Si no estás familiarizado con las matrices, piensa en ello como un sistema de ecuaciones lineales: estás buscando variables $x_1,...,x_{2000}$ de tal manera que la media de los primeros treinta produce $y_1$ la media de los segundos treinta da como resultado $y_2$ y así sucesivamente.

El problema del sistema de ecuaciones (y de la matriz) es que tiene más incógnitas que ecuaciones. Por lo tanto, no puede identificar de forma única sus incógnitas $x_1,...,x_n$ . La razón intuitiva es que se pierden dimensiones al promediar, porque las primeras treinta dimensiones de $\mathbf x$ no obtienen un elemento correspondiente en $\mathbf y$ ya que no se puede desplazar la ventana de promediación fuera de $\mathbf x$ .

Una forma de hacer $A$ o, lo que es lo mismo, el sistema de ecuaciones, solucionable es llegar a $30$ más ecuaciones (o $30$ más filas para $A$ ) que proporcionan información adicional (son linealmente independientes de todas las demás filas de $A$ ).

Otra forma, quizá más sencilla, es utilizar el pseudoinverso $A^\dagger$ de $A$ . Esto genera un vector $\mathbf z = A^\dagger\mathbf y$ que tiene la misma dimensión que $\mathbf x$ y que tiene la propiedad de minimizar la distancia cuadrática entre $\mathbf y$ y $A\mathbf z$ (ver wikipedia ).

Esto parece funcionar bastante bien. Aquí hay un ejemplo en el que dibujé $2000$ ejemplos a partir de una distribución gaussiana, añadimos cinco, los promediamos y reconstruimos el $\mathbf x$ a través del pseudoinverso.

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

Aquí estaría el código 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 eso ayude.

Answer 2

+1 a la respuesta de fabee, que es completa. Sólo un apunte para traducirlo a R, en base a los paquetes que he encontrado para hacer las operaciones en cuestión. En mi caso, tenía datos que son las previsiones de temperatura de la NOAA en base a tres meses: Ene-Feb-Mar, Feb-Mar-Abr, Mar-Abr-Mayo, etc, y quería descomponerlos en valores mensuales (aproximados), asumiendo que la temperatura de cada periodo de tres meses es esencialmente una media.

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 me funciona muy bien. Gracias @fabee.

EDIT: OK, volviendo a traducir mi R a Python, lo consigo:

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 de depuración que la versión de 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).

Answer 3

Esto está muy relacionado con esta pregunta cumsum con desplazamiento de n Pregunté en SO.

También respondí en SO la misma pregunta como este pero ha sido cerrado así que incluyo aquí la respuesta de nuevo porque creo que está más enfocado en la implementación del software que desde la comprensión matemática (aunque creo que son equivalentes matemáticamente).

El pregunta preguntó lo mismo, cómo invertir la media móvil, a.k.a en pandas como rolling mean.

La muestra de código del pregunta :

import numpy as np
import pandas as pd
import matplotlib.pylab as plt
np.random.seed(100)
data = np.random.rand(200,3)

df = pd.DataFrame(data)
df.columns = ['a', 'b', 'y']

df['y_roll'] = df['y'].rolling(10).mean()
df['y_roll_predicted'] = df['y_roll'].apply(lambda x: x + np.random.rand()/20)

Así que, cómo obtener df['y'] de vuelta de df['y_roll'] ? y aplicar el mismo método a df['y_roll_predicted']

Con esta función cumsum_shift(n) que tienes que pensar en ello como el inversa del método pandas/numpy diff(periods = n) , puede invertir la media móvil hasta la constante si no tiene los valores iniciales.

La definición de cumsum_shift(n) que generaliza el cumsum() que es éste con n = 1 (n se llama shift en el código):

def cumsum_shift(s, shift = 1, init_values = [0]):
    s_cumsum = pd.Series(np.zeros(len(s)))
    for i in range(shift):
        s_cumsum.iloc[i] = init_values[i]
    for i in range(shift,len(s)):
        s_cumsum.iloc[i] = s_cumsum.iloc[i-shift] + s.iloc[i]
    return s_cumsum

Entonces, suponiendo que el tamaño de la ventana es de 10 win_size = 10 entonces si se multiplica por 10 el diff 'ed de la media rodante y luego "cumsum shift it" con un desplazamiento de 10, se obtiene la serie original hasta los valores iniciales.

El código:

win_size = 10
s_diffed = win_size * df['y_roll'].diff()
df['y_unrolled'] = cumsum_shift(s=s_diffed, shift = win_size, init_values= df['y'].values[:win_size])

Este código recupera exactamente y de y_roll porque tienes los valores iniciales.

Puedes ver al trazarlo (en mi caso con plotly) que y y y_unrolled son exactamente iguales (sólo el rojo).

Ahora haciendo lo mismo a y_roll_predicted para obtener y_predicted_unrolled .

Código:

win_size = 10
s_diffed = win_size * df['y_roll_predicted'].diff()
df['y_predicted_unrolled'] = cumsum_shift(s=s_diffed, shift = win_size, init_values= df['y'].values[:win_size])

En este caso el resultado no es exactamente el mismo, fíjate en que el valores iniciales son de y y luego y_roll_predicted incorporar el ruido a y_roll por lo que el "desenrollado" no puede recuperar exactamente el original.

Aquí un gráfico ampliado en un rango más pequeño para verlo mejor:

Espero que esto pueda ayudar a alguien.

Answer 4

Gonzalo, Estoy usando tu función cumsum_shift en mi df grande (400.000 puntos) pero tengo problemas cuando cambio el win_size. La figura de abajo es para win_size=12,000 y puedo ver algunos picos al final de cada win_size. Para mi problema actual necesito usar win_size> 40.000. ¿Tiene alguna idea de la restricción de su función basada en el win_size? Gracias de antemano

Answer 5

La respuesta de Fabee fue completa. Sólo estoy añadiendo una función genérica que se puede utilizar en Python que he creado y probado para mis proyectos (con un código de ejemplo)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def reconstruct_orig(sm_x:np.ndarray, win_size:int=7):
    """reconstructing from original data

    Args:
        sm_x (np.ndarray): smoothed array (remove any NaN from the edge)
        win_size (int, optional): moving average window size. Defaults to 7.

    Returns:
        [type]: [description]
    """    '''
    '''
    arr_size = sm_x.shape[0]+win_size
    # get A and its inverse     
    A = (np.tril(np.ones((arr_size,arr_size)),-1) - np.tril(np.ones((arr_size,arr_size)),-(win_size+1)))/win_size
    A = A[win_size:,:]
    pA = np.linalg.pinv(A) #pseudo inverse
    return np.dot(pA, sm_x)

if __name__=="__main__":
    # np.random.seed(1)
    nmax= 100
    t=np.linspace(0,10,num=nmax)
    raw_x = pd.Series(np.sin(t)+ 0.2*np.random.normal(0,1, size=nmax)) # create original data
    sm_x = raw_x.rolling(7, center=False).mean().dropna() # smooth data
    re_x = reconstruct_orig(sm_x, win_size=7)          # reconstruct data

    plt.plot(raw_x,'x',label='original x')
    plt.plot(sm_x,label='averaged x')
    plt.plot(re_x,'.', label='reconstructed x')
    plt.legend()
    plt.show()