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?

  • 1
    $\begingroup$ 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? $\endgroup$ – whuber Aug 20 '13 at 20:17
  • $\begingroup$ @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! $\endgroup$ – user16679 Aug 20 '13 at 20:26
  • 1
    $\begingroup$ (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. $\endgroup$ – fabee Aug 21 '13 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.)

  • $\begingroup$ @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. $\endgroup$ – Wayne Aug 22 '13 at 19:30
  • $\begingroup$ OK, see if it works for you. If so, you could pick my answer ;-) $\endgroup$ – Wayne Aug 26 '13 at 2:02
  • $\begingroup$ 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... $\endgroup$ – wordsforthewise May 31 '17 at 2:13
  • $\begingroup$ Also, Jupyter notebooks for Python interactivity, or IPython. $\endgroup$ – wordsforthewise May 31 '17 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.

reconstruction of original signal from moving average using 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')

Hope that helps.

  • $\begingroup$ 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. $\endgroup$ – gdelfino Jul 26 '17 at 17:24
  • $\begingroup$ 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? $\endgroup$ – fabee Jul 28 '17 at 5:42
  • $\begingroup$ 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: $\endgroup$ – gdelfino Jul 28 '17 at 14:35
  • $\begingroup$ {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} $\endgroup$ – gdelfino Jul 28 '17 at 14:35
  • $\begingroup$ 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. $\endgroup$ – gdelfino Jul 31 '17 at 13:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy