# Extract data points from moving average?

Is it possible to extract data points from moving average data?

In other words, if a set of data only has simple moving averages of the previous 30 points, is it possible to extract the original data points?

If so, how?

• The answer is a qualified yes, but the exact procedure depends on how the initial segment of data is treated. If it is simply dropped, then you have effectively lost 15 pieces of data, leaving you with an undeterdetermined system of linear equations. The upshot is that there exist many valid answers in general, but you can still make some progress if either (a) shorter windows (or some such procedure) are used for the initial 15 moving averages or (b) you can specify additional constraints on the solution (about 15 dimensions' worth of constraints...). What situation are you in?
– whuber
Aug 20, 2013 at 20:17
• @whuber Thank you very much for looking! I have 2,000 points. The first MA point is most likely an average of the first 30 original points. Accuracy is second to a generally correct result, most specifically good guesses at the most "recent" points. Can you recommend a relatively simple method? Thanks in advance!
– user16679
Aug 20, 2013 at 20:26
• (if you take more than five minutes to write a comment ...). What I wanted to write is that you can think of the averaging as a matrix multiplication. The rows in the middle will have 1/30*[1 1 1 ...] before the diagonal. The question is, how do you deal with points at the borders of your vector in order to make the matrix invertible. You can do this by assuming that they are a result of averaging over less elements or you think about other constraints. Note that while a matrix inversion is an easy way to understand it, it is not the most efficient. You probably want to use an FFT to do that. Aug 21, 2013 at 6:06

+1 to fabee's answer, which is complete. Just a note to translate it into R, based on the packages that I've found to do the operations at hand. In my case, I had data that is NOAA temperature forecasts on a three-month basis: Jan-Feb-Mar, Feb-Mar-Apr, Mar-Apr-May, etc, and I wanted to break it out into (approximate) monthly values, assuming that each three-month period's temperature is essentially an average.

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


Which works great for me. Thanks @fabee.

EDIT: OK, back-translating my R to Python, I get:

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


(Which took a lot longer to debug than the R version. First because I'm not as familiar with Python as with R, but also because R is much more usable interactively.)

• @Gracchus: Sorry, not a C++ guy, but you may find what you need in the Armadillo C++ linear algebra library (arma.sourceforge.net), which is also available in R via the RcppArmadillo package. Aug 22, 2013 at 19:30
• OK, see if it works for you. If so, you could pick my answer ;-) Aug 26, 2013 at 2:02
• FYI best practices in Python are to do absolute imports: python.org/dev/peps/pep-0008/#imports which makes it so much easier to read other people's code, because you actually know where the functions come from instead of having to look up each one you don't know. Wish it was standard in R to do the same. Having to lookup every little functions in someone else's code really grinds my gears... May 31, 2017 at 2:13
• Also, Jupyter notebooks for Python interactivity, or IPython. May 31, 2017 at 2:19

I try to put what whuber said into an answer. Let's say you have a large vector $\mathbf x$ with $n=2000$ entries. If you compute a moving average with a window of length $\ell=30$, you can write this as a vector matrix multiplication $\mathbf y = A\mathbf x$ of the vector $\mathbf x$ with the matrix

$$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)$$

which has $30$ ones which are shifted through as you advance through the rows until the $30$ ones hit the end of the matrix. Here the averaged vector $\mathbf y$ has 1970 dimensions. The matrix has $1970$ rows and $2000$ columns. Therefore, it is not invertible.

If you are not familiar with matrices, think about it as a linear equation system: you are searching for variables $x_1,...,x_{2000}$ such that the average over the first thirty yields $y_1$, the average over the second thirty yields $y_2$ and so on.

The problem with the equation system (and the matrix) is that it has more unknowns than equations. Therefore, you cannot uniquely identify your unknowns $x_1,...,x_n$. The intuitive reason is that you loose dimensions while averaging, because the first thirty dimensions of $\mathbf x$ don't get a corresponding element in $\mathbf y$ since you cannot shift the averaging window outside of $\mathbf x$.

One way to make $A$ or, equivalently the equation system, solvable is to come up with $30$ more equations (or $30$ more rows for $A$) that provide additional information (are linearly independent to all other rows of $A$).

Another, maybe easier, way is to use the pseudoinverse $A^\dagger$ of $A$. This generates a vector $\mathbf z = A^\dagger\mathbf y$ which has the same dimension as $\mathbf x$ and which has the property that it minimizes the quadratic distance between $\mathbf y$ and $A\mathbf z$ (see wikipedia).

This seems to work quite well. Here is an example where I drew $2000$ examples from a Gaussian distribution, added five, averaged them, and reconstructed the $\mathbf x$ via the pseudoinverse.

Many numerical programs offer pseudo-inverses (e.g. Matlab, numpy in python, etc.).

Here would be the python code to generate the signals from my example:

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()


Hope that helps.

• This is a great answer, but I think you are mistaken when you said that "it minimizes the quadratic distance between y and A.z". In fact y and A.z are the same thing. What is minimized is the norm of z which works well for the real world signals that I have tried, but is not so good if your original signal has many outliers. Jul 26, 2017 at 17:24
• I am not sure whether I follow. y and A.x are the same thing, but not y and A.z. It's true that it also minimizes the norm of z. I also don't see why it doesn't work for my examples. The blue and the red line match up pretty nicely. Am I missing something in your comment? Jul 28, 2017 at 5:42
• y is the moving average calculated from the original signal x by multiplying by A. This procedure give us a signal z which has the same moving average y. Therefore y=A.z. So only the norm of z gets minimized. If the original signal happen to have a large norm value, then the procedure will not give good results. An example signal with large norm value is below: Jul 28, 2017 at 14:35
• {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.8894, 20.4211, 43.4591, -33.2786, 13.5468, -45.2374, 10.3787, 35.8235, -41.5161, 19.9717, 43.0658, -33.7125, 13.0321} Jul 28, 2017 at 14:35
• Please use a windows size of 8 for the signal above. This way the filtered signal is very different in shape from the original signal. Jul 31, 2017 at 13:58

This is very related with this question cumsum with shift of n I asked in SO.

I also answered in SO the same question as this one but it has been closed so I include here the answer again because I think is more focus in the software implementation than from the mathematical understanding (even though I think they are equivalent mathematically).

The question asked the same thing, how to reverse the moving average, a.k.a in pandas as rolling mean.

The code sample of the question:

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)


So, how to obtain df['y'] back from df['y_roll']? and apply the same method to df['y_roll_predicted']

With this function cumsum_shift(n) which you have to think of it as the inverse of the pandas/numpy method diff(periods = n), you can reverse the moving average up to constant if you don't have the initial values.

The definition of cumsum_shift(n) that generalizes the cumsum() which is this one with n = 1 (n is called shift in the code):

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


Then assuming the size of the window is 10 win_size = 10 then if you multiply by 10 the diff'ed of the rolling mean and then "cumsum shift it" with a shift of 10, you obtain the original serie up to the intial values.

The code:

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])


This code recovers exactly y from y_roll because you have the initial values.

You can see it plotting it (in my case with plotly) that y and y_unrolled are exactly the same (just the red one).

Now doing the same thing to y_roll_predicted to obtain y_predicted_unrolled.

Code:

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])


In this case the result are not exactly the same, notice how the initial values are from y and then y_roll_predicted incorporate noise to y_roll so the "unrolling" cannot recover exactly the original one.

Here a plot zoomed in in a smaller range to see it better:

Hope this can help somebody.

Gonzalo, I'm using your cumsum_shift function in my large df (400,000 points) but I have problems when I change the win_size. Figure below is for win_size=12,000 and I can see some spikes at the end of each win_size. For my current problem I need to use win_size> 40,000. Do you have any idea of restriction of your function based on the win_size? Thanks in advance

• If you have a new question, please ask it by clicking the Ask Question button. Include a link to this question if it helps provide context. - From Review Oct 17, 2021 at 18:17

fabee's answer was complete. I am just adding a generic function that can be used in Python that I've created and tested for my projects (with a sample code)

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()