RE: vectorized moving average?

[Tim Rueth]

Well, I can't figure out how you came up with the initial condition:  (data(1)*(1-alpha)), but it agrees perfectly with the for-loop version shown below.  Thanks!

 

--Tim

---

From: address@hidden [mailto:address@hidden On Behalf Of James Sherman Jr.
Sent: Friday, May 14, 2010 6:57 AM
To: address@hidden
Cc: Octave-ML
Subject: Re: vectorized moving average?

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