# 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 '13 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 '13 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. –  fabee Aug 21 '13 at 6:06

## 2 Answers

+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 then 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. –  Wayne Aug 22 '13 at 19:30
OK, see if it works for you. If so, you could pick my answer ;-) –  Wayne Aug 26 '13 at 2:02
thank you so much for your hard work! –  user16679 Aug 26 '13 at 2:25

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)

# 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.

-