Re: vectorized moving average?

[James Sherman Jr.]

So, this bugged me, so I looked a bit at the filter function, and I think I found where the mistake was in my initial suggestion.  The initial condition vector has to do with the internal states of the filter not the negative time outputs of the filter (at least not directly), so to get what I think is exactly what your code with the for loop, the filter line should be:

avg = filter(alpha, [1 alpha-1], data, data(1)*(1-alpha));

It is rather unintuitive why the 1-alpha term needs to be there, and I don't know if there's much interest in it, but it shouldn't be that hard (probably I just need to crack open my signals and systems book) to write a function to calculate the those initial conditions that the filter function expects just giving the outputs and inputs from negative time.

On Thu, May 13, 2010 at 8:38 PM, Tim Rueth <address@hidden> wrote:

The last instruction with "long_ma" should have read: "avg = avg(n+1 :
end);" which effectively trims off the computed values from negative time.
But, as you say, it looks like I didn't need to do this because the history
is completely captured in avg(1) = data(1), so no need to compute a "run-in"
time.  Thanks Francesco.

Sherman had found that I can set the initial condition by specifying a 4th
parameter in filter() equal to the first data point.  I tried this, and got
very similar (but not quite exact) results when compared to the for-loop
below with no negative time values.  But this small difference dissipated
within "ndays" and isn't a big deal.  Thanks Sherman.

In summary, to calculate the exponential moving average of "data" for
"ndays", the following code:

 alpha = 2/(ndays+1);
 n = length(data);
 avg = zeros(n,1);
 avg(1) = data(1);

 for i = 2 : n
    ao = avg(i-1);
    avg(i) = ao + alpha*(data(i) - ao);

 endfor;

...is close, but not quite equal to:

 alpha = 2/(ndays+1);
 avg = filter(alpha, [1 alpha-1], data, data(1));

...for roughly the first ndays of avg.

--Tim

> -----Original Message-----
> From: Francesco Potortì [mailto:address@hidden]
> Sent: Wednesday, May 12, 2010 11:22 PM
> To: address@hidden

> Cc: 'Octave-ML'; 'James Sherman Jr.'
> Subject: Re: vectorized moving average?
>

> >Your filter code below works just fine when compared to what
> I had been
> >doing, except for a number of initial days, due to what values are
> >assumed in negative time.  I had been using the following code:
> >
> ># "ndays" is the number of days to be used when computing the
> >exponential moving average of "data" (data is a column vector)
> > data = "" ndays, 1); data];   # repeat
> data(1) ndays times at
> >the beginning of data for negative time values  alpha =
> 2/(ndays+1);  n
> >= length(data);  avg = zeros(n,1);
> > avg(1) = data(1);
>
> The above instruction is all you need to "invent" past memory
> for negative values.  You should do the same for the filter
> function, but I could not say how to do it offhand.
>
> > for i = 2 : n
> >     ao = avg(i-1);
> >     avg(i) = ao + alpha*(data(i) - ao); endfor
> >
> ># trim off run-in period for negative time values
> >   long_ma = long_ma(lma_days+1 : end);
>
> I don't understand the above instruction.  What is long_ma?
>
> >For small values of ndays, the number of initial days where
> there's a
> >discrepancy with your filter() implementation is minimal, but for
> >larger values of ndays, the number of initial days of
> discrepancy grows
> >(obviously, due to the nature of an exponential MA having a
> long-tail
> >memory).  Note, I add similar negative time values to the
> front of the
> >vector when using
> >filter() as well.  I'm just not sure what is the convention when it
> >comes to calculating exponential moving averages for points
> in "data" where "ndays"
> >reaches back into negative time.  Thanks again.
>
> --
> Francesco Potortì (ricercatore)        Voice: +39 050 315
> 3058 (op.2111)
> ISTI - Area della ricerca CNR          Fax:   +39 050 315 2040
> via G. Moruzzi 1, I-56124 Pisa         Email: address@hidden
> (entrance 20, 1st floor, room C71)     Web:   http://fly.isti.cnr.it/
>