Does anyone know how
to take an n-day weighted moving average of a vector without using a
for-loop? I looked at the M code for movavg() and it uses a for-loop, so
I'm guessing there probably isn't a way, but I thought I'd check.
Thanks.
--Tim
_______________________________________________ Help-octave mailing list [hidden email] https://www-old.cae.wisc.edu/mailman/listinfo/help-octave |
conv is also an m file, but it only has a few ifs in. then it calls
filter to get the job done. which is an oct file. Andy On Thu, May 6, 2010 at 6:28 AM, Tim Rueth <[hidden email]> wrote: > Does anyone know how to take an n-day weighted moving average of a vector > without using a for-loop? I looked at the M code for movavg() and it uses a > for-loop, so I'm guessing there probably isn't a way, but I thought I'd > check. Thanks. > > --Tim > _______________________________________________ > Help-octave mailing list > [hidden email] > https://www-old.cae.wisc.edu/mailman/listinfo/help-octave > > -- /* andy buckle */ _______________________________________________ Help-octave mailing list [hidden email] https://www-old.cae.wisc.edu/mailman/listinfo/help-octave |
In reply to this post by Tim R
>Does anyone know how to take an n-day weighted moving average of a vector
>without using a for-loop? Use the filter function. > I looked at the M code for movavg() and it uses a >for-loop, so I'm guessing there probably isn't a way, but I thought I'd >check. Thanks. I don't know why it does so, but it looks very inefficient, at least for long vectors. -- 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: [hidden email] (entrance 20, 1st floor, room C71) Web: http://fly.isti.cnr.it/ _______________________________________________ Help-octave mailing list [hidden email] https://www-old.cae.wisc.edu/mailman/listinfo/help-octave |
On Thu, May 6, 2010 at 11:40 AM, Francesco Potortì <[hidden email]> wrote:
>Does anyone know how to take an n-day weighted moving average of a vector Or fftconv /Fredrik _______________________________________________ Help-octave mailing list [hidden email] https://www-old.cae.wisc.edu/mailman/listinfo/help-octave |
On Thu, May 6, 2010 at 12:45 PM, Fredrik Lingvall <[hidden email]> wrote:
> > > On Thu, May 6, 2010 at 11:40 AM, Francesco Potortì <[hidden email]> > wrote: >> >> >Does anyone know how to take an n-day weighted moving average of a vector >> >without using a for-loop? >> >> Use the filter function. > > Or fftconv > Or, for unweighted moving average, probably the fastest way is to compute differences of a cumulative sum: y = cumsum (x); n = rows (x); avg = (y(k+1:n) - y(1:n-k+1)) / k; unfortunately, this way is also probably the most inaccurate one :) -- RNDr. Jaroslav Hajek, PhD computing expert & GNU Octave developer Aeronautical Research and Test Institute (VZLU) Prague, Czech Republic url: www.highegg.matfyz.cz _______________________________________________ Help-octave mailing list [hidden email] https://www-old.cae.wisc.edu/mailman/listinfo/help-octave |
In reply to this post by andy buckle
I looked at both conv() and filter(), but can't figure out how to do a
moving average with them. Perhaps I'm not understanding the functions of the input vars correctly. Let's say I have an array, a = rand(1,100). Can you tell me how I'd use conv() and filter() to take, say the 10-day moving average of it, with a weighting of 0.5? Thanks! > -----Original Message----- > From: Andy Buckle [mailto:[hidden email]] > Sent: Thursday, May 06, 2010 12:06 AM > To: [hidden email] > Cc: [hidden email] > Subject: Re: vectorized moving average? > > conv is also an m file, but it only has a few ifs in. then it > calls filter to get the job done. which is an oct file. > > Andy > > On Thu, May 6, 2010 at 6:28 AM, Tim Rueth <[hidden email]> wrote: > > Does anyone know how to take an n-day weighted moving average of a > > vector without using a for-loop? I looked at the M code > for movavg() > > and it uses a for-loop, so I'm guessing there probably isn't a way, > > but I thought I'd check. Thanks. > > > > --Tim > > _______________________________________________ > > Help-octave mailing list > > [hidden email] > > https://www-old.cae.wisc.edu/mailman/listinfo/help-octave > > > > > > > > -- > /* andy buckle */ > _______________________________________________ Help-octave mailing list [hidden email] https://www-old.cae.wisc.edu/mailman/listinfo/help-octave |
I'm not sure what you mean by weighting of 0.5, but to do a simple 10-day average, it'd be
x = rand(1, 100); y = filter(1/10*ones(1, 10), 1, x) This assumes that the values at negative time (x(0), x(-1), etc) are all zero. So, for example, the first value of y would be x(1)/10.
On Fri, May 7, 2010 at 3:33 PM, Tim Rueth <[hidden email]> wrote: I looked at both conv() and filter(), but can't figure out how to do a _______________________________________________ Help-octave mailing list [hidden email] https://www-old.cae.wisc.edu/mailman/listinfo/help-octave |
Thanks for showing how to use filter() to do a simple moving
average. I implemented your code, and it agrees with movavg(x,10,10,0),
which calculates a 10-day simple moving average. There's only a difference
in the first 9 numbers due to assumed values of negative time (movavg
calculates a "run-in" period). As you might recall, I'm trying to use
filter() to avoid movavg()'s for-loop.
Now, what I'm trying to do is a weighted moving average,
identical to the "alpha" parameter of movavg(). When alpha=0, it's a
simple moving average, and agrees with filter(). If I change alpha to 1,
I'm suppose to get a linear MA. Here's the code in movavg.m that does the
weighting (lead is the number of days to average, equal to 10 in the above
case):
lead = (1:lead).^alpha;
## Adjust the weights to equal 1 lead = lead / sum (lead); So, for a 10-day linearly-weighted moving
average (lead = 10, alpha = 1),
the previous 9 days and current day
should be weighted as follows: 1/55, 2/55, 3/55...10/55, with the greatest weight (10/55) applied to the current day. So, I tried a simple test case with just a 2-day MA on
a 6-element vector:
ma_days = 2;
alpha = 1; len = 6; a = rand(1,len); # Calculate MA using
movavg()
ma = movavg(a,ma_days,ma_days,alpha); # Calculate MA using
filter()
sweep = (1:ma_days).^alpha; norm_sweep = sweep/sum(sweep); f = filter(norm_sweep,1,a); disp(ma);
disp(f); The movavg() and filter() results are similar, but not
equal. I'm guessing I don't have the arguments for filter() correct, but I
can't figure out what I did wrong. Particularly, I'm not sure what the
second argument of filter() is supposed to do.
Help?
_______________________________________________ Help-octave mailing list [hidden email] https://www-old.cae.wisc.edu/mailman/listinfo/help-octave |
I can't double-check right now, but I think it is that filter is setup
so that it does (in this simplified case), let b be the coefficients and has length m, y(n) = b(1)*x(n) + b(2)*x(n-1) ... + b(m)*x(n-m+1) So, the first value of your coefficient vector (sweep as you call it) should be the value you want to weigh the "current" value of x (what you call a). Thus, long story short, I believe that your sweep value should be: sweep = (ma_days:-1:1).^alpha; norm_sweep = sweep/sum(sweep); Just reversed of what you put in. Also, as a side note, the reason why there wasn't much difference between movavg and filter for what you put in, is because you just gave it uniform noise, so any average function is really just approximating the average of the noise. On Tue, May 11, 2010 at 1:41 PM, Tim Rueth <[hidden email]> wrote: > > Thanks for showing how to use filter() to do a simple moving average. I implemented your code, and it agrees with movavg(x,10,10,0), which calculates a 10-day simple moving average. There's only a difference in the first 9 numbers due to assumed values of negative time (movavg calculates a "run-in" period). As you might recall, I'm trying to use filter() to avoid movavg()'s for-loop. > > Now, what I'm trying to do is a weighted moving average, identical to the "alpha" parameter of movavg(). When alpha=0, it's a simple moving average, and agrees with filter(). If I change alpha to 1, I'm suppose to get a linear MA. Here's the code in movavg.m that does the weighting (lead is the number of days to average, equal to 10 in the above case): > > lead = (1:lead).^alpha; > ## Adjust the weights to equal 1 > lead = lead / sum (lead); > > So, for a 10-day linearly-weighted moving average (lead = 10, alpha = 1), the previous 9 days and current day should be weighted as follows: 1/55, 2/55, 3/55...10/55, with the greatest weight (10/55) applied to the current day. So, I tried a simple test case with just a 2-day MA on a 6-element vector: > > ma_days = 2; > alpha = 1; > len = 6; > a = rand(1,len); > > # Calculate MA using movavg() > ma = movavg(a,ma_days,ma_days,alpha); > > # Calculate MA using filter() > sweep = (1:ma_days).^alpha; > norm_sweep = sweep/sum(sweep); > f = filter(norm_sweep,1,a); > > disp(ma); > disp(f); > > The movavg() and filter() results are similar, but not equal. I'm guessing I don't have the arguments for filter() correct, but I can't figure out what I did wrong. Particularly, I'm not sure what the second argument of filter() is supposed to do. Help? > > ________________________________ > From: [hidden email] [mailto:[hidden email]] On Behalf Of James Sherman Jr. > Sent: Friday, May 07, 2010 12:58 PM > To: [hidden email] > Cc: Octave-ML > Subject: Re: vectorized moving average? > > I'm not sure what you mean by weighting of 0.5, but to do a simple 10-day average, it'd be > x = rand(1, 100); > y = filter(1/10*ones(1, 10), 1, x) > This assumes that the values at negative time (x(0), x(-1), etc) are all zero. So, for example, the first value of y would be x(1)/10. > On Fri, May 7, 2010 at 3:33 PM, Tim Rueth <[hidden email]> wrote: >> >> I looked at both conv() and filter(), but can't figure out how to do a >> moving average with them. Perhaps I'm not understanding the functions of >> the input vars correctly. >> >> Let's say I have an array, a = rand(1,100). Can you tell me how I'd use >> conv() and filter() to take, say the 10-day moving average of it, with a >> weighting of 0.5? >> >> Thanks! >> >> >> > -----Original Message----- >> > From: Andy Buckle [mailto:[hidden email]] >> > Sent: Thursday, May 06, 2010 12:06 AM >> > To: [hidden email] >> > Cc: [hidden email] >> > Subject: Re: vectorized moving average? >> > >> > conv is also an m file, but it only has a few ifs in. then it >> > calls filter to get the job done. which is an oct file. >> > >> > Andy >> > >> > On Thu, May 6, 2010 at 6:28 AM, Tim Rueth <[hidden email]> wrote: >> > > Does anyone know how to take an n-day weighted moving average of a >> > > vector without using a for-loop? I looked at the M code >> > for movavg() >> > > and it uses a for-loop, so I'm guessing there probably isn't a way, >> > > but I thought I'd check. Thanks. >> > > >> > > --Tim >> > > _______________________________________________ >> > > Help-octave mailing list >> > > [hidden email] >> > > https://www-old.cae.wisc.edu/mailman/listinfo/help-octave >> > > >> > > >> > >> > >> > >> > -- >> > /* andy buckle */ >> > >> >> >> _______________________________________________ >> Help-octave mailing list >> [hidden email] >> https://www-old.cae.wisc.edu/mailman/listinfo/help-octave > _______________________________________________ Help-octave mailing list [hidden email] https://www-old.cae.wisc.edu/mailman/listinfo/help-octave |
Yep, worked great, thanks!
Now, I just have to figure out how to do an exponential MA using filter(). (I think movavg() has a bug in it for the exponential case, where alpha = 'e'.) --Tim > -----Original Message----- > From: [hidden email] [mailto:[hidden email]] On > Behalf Of James Sherman Jr. > Sent: Tuesday, May 11, 2010 11:01 AM > To: [hidden email] > Cc: Octave-ML > Subject: Re: vectorized moving average? > > I can't double-check right now, but I think it is that filter > is setup so that it does (in this simplified case), let b be > the coefficients and has length m, > > y(n) = b(1)*x(n) + b(2)*x(n-1) ... + b(m)*x(n-m+1) > > So, the first value of your coefficient vector (sweep as you > call it) should be the value you want to weigh the "current" > value of x (what you call a). Thus, long story short, I > believe that your sweep value should be: > > sweep = (ma_days:-1:1).^alpha; > norm_sweep = sweep/sum(sweep); > > Just reversed of what you put in. Also, as a side note, the > reason why there wasn't much difference between movavg and > filter for what you put in, is because you just gave it > uniform noise, so any average function is really just > approximating the average of the noise. > > On Tue, May 11, 2010 at 1:41 PM, Tim Rueth <[hidden email]> wrote: > > > > Thanks for showing how to use filter() to do a simple > moving average. I implemented your code, and it agrees with > movavg(x,10,10,0), which calculates a 10-day simple moving > average. There's only a difference in the first 9 numbers > due to assumed values of negative time (movavg calculates a > "run-in" period). As you might recall, I'm trying to use > filter() to avoid movavg()'s for-loop. > > > > Now, what I'm trying to do is a weighted moving average, > identical to the "alpha" parameter of movavg(). When > alpha=0, it's a simple moving average, and agrees with > filter(). If I change alpha to 1, I'm suppose to get a > linear MA. Here's the code in movavg.m that does the > weighting (lead is the number of days to average, equal to 10 > in the above case): > > > > lead = (1:lead).^alpha; > > ## Adjust the weights to equal 1 > > lead = lead / sum (lead); > > > > So, for a 10-day linearly-weighted moving average (lead = > 10, alpha = 1), the previous 9 days and current day should be > weighted as follows: 1/55, 2/55, 3/55...10/55, with the > greatest weight (10/55) applied to the current day. So, I > tried a simple test case with just a 2-day MA on a 6-element vector: > > > > ma_days = 2; > > alpha = 1; > > len = 6; > > a = rand(1,len); > > > > # Calculate MA using movavg() > > ma = movavg(a,ma_days,ma_days,alpha); > > > > # Calculate MA using filter() > > sweep = (1:ma_days).^alpha; > > norm_sweep = sweep/sum(sweep); > > f = filter(norm_sweep,1,a); > > > > disp(ma); > > disp(f); > > > > The movavg() and filter() results are similar, but not > equal. I'm guessing I don't have the arguments for filter() > correct, but I can't figure out what I did wrong. > Particularly, I'm not sure what the second argument of > filter() is supposed to do. Help? > > > > ________________________________ > > From: [hidden email] [mailto:[hidden email]] On > Behalf Of James Sherman Jr. > > Sent: Friday, May 07, 2010 12:58 PM > > To: [hidden email] > > Cc: Octave-ML > > Subject: Re: vectorized moving average? > > > > I'm not sure what you mean by weighting of 0.5, but to do a simple > > 10-day average, it'd be x = rand(1, 100); y = > filter(1/10*ones(1, 10), > > 1, x) This assumes that the values at negative time (x(0), > x(-1), etc) > > are all zero. So, for example, the first value of y would > be x(1)/10. > > On Fri, May 7, 2010 at 3:33 PM, Tim Rueth > <[hidden email]> wrote: > >> > >> I looked at both conv() and filter(), but can't figure out > how to do > >> a moving average with them. Perhaps I'm not understanding the > >> functions of the input vars correctly. > >> > >> Let's say I have an array, a = rand(1,100). Can you tell > me how I'd > >> use > >> conv() and filter() to take, say the 10-day moving average of it, > >> with a weighting of 0.5? > >> > >> Thanks! > >> > >> > >> > -----Original Message----- > >> > From: Andy Buckle [mailto:[hidden email]] > >> > Sent: Thursday, May 06, 2010 12:06 AM > >> > To: [hidden email] > >> > Cc: [hidden email] > >> > Subject: Re: vectorized moving average? > >> > > >> > conv is also an m file, but it only has a few ifs in. > then it calls > >> > filter to get the job done. which is an oct file. > >> > > >> > Andy > >> > > >> > On Thu, May 6, 2010 at 6:28 AM, Tim Rueth > <[hidden email]> wrote: > >> > > Does anyone know how to take an n-day weighted moving > average of > >> > > a vector without using a for-loop? I looked at the M code > >> > for movavg() > >> > > and it uses a for-loop, so I'm guessing there probably isn't a > >> > > way, but I thought I'd check. Thanks. > >> > > > >> > > --Tim > >> > > _______________________________________________ > >> > > Help-octave mailing list > >> > > [hidden email] > >> > > https://www-old.cae.wisc.edu/mailman/listinfo/help-octave > >> > > > >> > > > >> > > >> > > >> > > >> > -- > >> > /* andy buckle */ > >> > > >> > >> > >> _______________________________________________ > >> Help-octave mailing list > >> [hidden email] > >> https://www-old.cae.wisc.edu/mailman/listinfo/help-octave > > > _______________________________________________ Help-octave mailing list [hidden email] https://www-old.cae.wisc.edu/mailman/listinfo/help-octave |
Well, if you want to use the formula found in Wikipedia exponential moving average, its not as you explain it, because an exponential moving average is fundamentally different than n-day weighted average because it has an infinite "memory" or from my EE background, it has an infinite impulse response. (I don't know what your background is).
But don't fear! The filter command can handle this as well (quite easily in fact). You just need to use the equation found in in the wikipedia article that shows (let s_t be the sum at time t, and x_t be the value at time t), then s_t = s_(t-1) + alpha*(x_t - s_(t-1)) or rearranged, it looks like s_t = alpha*x_t - (alpha - 1)*s_(t-1) so comparing those values to what filter expects, we get: b = alpha; a = [1, alpha-1]; s = filter(b, a, x); I'm not 100% I haven't made a mistake here, but this is the general gist of things. Hope this helps. On Tue, May 11, 2010 at 4:28 PM, Tim Rueth <[hidden email]> wrote: > Yep, worked great, thanks! > > Now, I just have to figure out how to do an exponential MA using filter(). > (I think movavg() has a bug in it for the exponential case, where alpha = > 'e'.) > > --Tim > > >> -----Original Message----- >> From: [hidden email] [mailto:[hidden email]] On >> Behalf Of James Sherman Jr. >> Sent: Tuesday, May 11, 2010 11:01 AM >> To: [hidden email] >> Cc: Octave-ML >> Subject: Re: vectorized moving average? >> >> I can't double-check right now, but I think it is that filter >> is setup so that it does (in this simplified case), let b be >> the coefficients and has length m, >> >> y(n) = b(1)*x(n) + b(2)*x(n-1) ... + b(m)*x(n-m+1) >> >> So, the first value of your coefficient vector (sweep as you >> call it) should be the value you want to weigh the "current" >> value of x (what you call a). Thus, long story short, I >> believe that your sweep value should be: >> >> sweep = (ma_days:-1:1).^alpha; >> norm_sweep = sweep/sum(sweep); >> >> Just reversed of what you put in. Also, as a side note, the >> reason why there wasn't much difference between movavg and >> filter for what you put in, is because you just gave it >> uniform noise, so any average function is really just >> approximating the average of the noise. >> >> On Tue, May 11, 2010 at 1:41 PM, Tim Rueth <[hidden email]> wrote: >> > >> > Thanks for showing how to use filter() to do a simple >> moving average. I implemented your code, and it agrees with >> movavg(x,10,10,0), which calculates a 10-day simple moving >> average. There's only a difference in the first 9 numbers >> due to assumed values of negative time (movavg calculates a >> "run-in" period). As you might recall, I'm trying to use >> filter() to avoid movavg()'s for-loop. >> > >> > Now, what I'm trying to do is a weighted moving average, >> identical to the "alpha" parameter of movavg(). When >> alpha=0, it's a simple moving average, and agrees with >> filter(). If I change alpha to 1, I'm suppose to get a >> linear MA. Here's the code in movavg.m that does the >> weighting (lead is the number of days to average, equal to 10 >> in the above case): >> > >> > lead = (1:lead).^alpha; >> > ## Adjust the weights to equal 1 >> > lead = lead / sum (lead); >> > >> > So, for a 10-day linearly-weighted moving average (lead = >> 10, alpha = 1), the previous 9 days and current day should be >> weighted as follows: 1/55, 2/55, 3/55...10/55, with the >> greatest weight (10/55) applied to the current day. So, I >> tried a simple test case with just a 2-day MA on a 6-element vector: >> > >> > ma_days = 2; >> > alpha = 1; >> > len = 6; >> > a = rand(1,len); >> > >> > # Calculate MA using movavg() >> > ma = movavg(a,ma_days,ma_days,alpha); >> > >> > # Calculate MA using filter() >> > sweep = (1:ma_days).^alpha; >> > norm_sweep = sweep/sum(sweep); >> > f = filter(norm_sweep,1,a); >> > >> > disp(ma); >> > disp(f); >> > >> > The movavg() and filter() results are similar, but not >> equal. I'm guessing I don't have the arguments for filter() >> correct, but I can't figure out what I did wrong. >> Particularly, I'm not sure what the second argument of >> filter() is supposed to do. Help? >> > >> > ________________________________ >> > From: [hidden email] [mailto:[hidden email]] On >> Behalf Of James Sherman Jr. >> > Sent: Friday, May 07, 2010 12:58 PM >> > To: [hidden email] >> > Cc: Octave-ML >> > Subject: Re: vectorized moving average? >> > >> > I'm not sure what you mean by weighting of 0.5, but to do a simple >> > 10-day average, it'd be x = rand(1, 100); y = >> filter(1/10*ones(1, 10), >> > 1, x) This assumes that the values at negative time (x(0), >> x(-1), etc) >> > are all zero. So, for example, the first value of y would >> be x(1)/10. >> > On Fri, May 7, 2010 at 3:33 PM, Tim Rueth >> <[hidden email]> wrote: >> >> >> >> I looked at both conv() and filter(), but can't figure out >> how to do >> >> a moving average with them. Perhaps I'm not understanding the >> >> functions of the input vars correctly. >> >> >> >> Let's say I have an array, a = rand(1,100). Can you tell >> me how I'd >> >> use >> >> conv() and filter() to take, say the 10-day moving average of it, >> >> with a weighting of 0.5? >> >> >> >> Thanks! >> >> >> >> >> >> > -----Original Message----- >> >> > From: Andy Buckle [mailto:[hidden email]] >> >> > Sent: Thursday, May 06, 2010 12:06 AM >> >> > To: [hidden email] >> >> > Cc: [hidden email] >> >> > Subject: Re: vectorized moving average? >> >> > >> >> > conv is also an m file, but it only has a few ifs in. >> then it calls >> >> > filter to get the job done. which is an oct file. >> >> > >> >> > Andy >> >> > >> >> > On Thu, May 6, 2010 at 6:28 AM, Tim Rueth >> <[hidden email]> wrote: >> >> > > Does anyone know how to take an n-day weighted moving >> average of >> >> > > a vector without using a for-loop? I looked at the M code >> >> > for movavg() >> >> > > and it uses a for-loop, so I'm guessing there probably isn't a >> >> > > way, but I thought I'd check. Thanks. >> >> > > >> >> > > --Tim >> >> > > _______________________________________________ >> >> > > Help-octave mailing list >> >> > > [hidden email] >> >> > > https://www-old.cae.wisc.edu/mailman/listinfo/help-octave >> >> > > >> >> > > >> >> > >> >> > >> >> > >> >> > -- >> >> > /* andy buckle */ >> >> > >> >> >> >> >> >> _______________________________________________ >> >> Help-octave mailing list >> >> [hidden email] >> >> https://www-old.cae.wisc.edu/mailman/listinfo/help-octave >> > >> > > _______________________________________________ Help-octave mailing list [hidden email] https://www-old.cae.wisc.edu/mailman/listinfo/help-octave |
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 = [repmat(data(1), 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); 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);
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.
_______________________________________________ Help-octave mailing list [hidden email] https://www-old.cae.wisc.edu/mailman/listinfo/help-octave |
>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 = [repmat(data(1), 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: [hidden email] (entrance 20, 1st floor, room C71) Web: http://fly.isti.cnr.it/ _______________________________________________ Help-octave mailing list [hidden email] https://www-old.cae.wisc.edu/mailman/listinfo/help-octave |
I can't check it currently, but if I remember correctly, the 4th argument to filter is initial conditions. So, something like if you want your initial condition to be the first value of data, I think the command would be:
b = alpha; a = [1, alpha-1]; s = filter(b, a, x, x(1)); It only needs to be one element in this case because the only initial condition you need is s_0. On Thu, May 13, 2010 at 2:21 AM, Francesco Potortì <[hidden email]> wrote:
_______________________________________________ Help-octave mailing list [hidden email] https://www-old.cae.wisc.edu/mailman/listinfo/help-octave |
In reply to this post by Francesco Potortì
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:[hidden email]] > Sent: Wednesday, May 12, 2010 11:22 PM > To: [hidden email] > 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 = [repmat(data(1), 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: [hidden email] > (entrance 20, 1st floor, room C71) Web: http://fly.isti.cnr.it/ > _______________________________________________ Help-octave mailing list [hidden email] https://www-old.cae.wisc.edu/mailman/listinfo/help-octave |
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 <[hidden email]> wrote: The last instruction with "long_ma" should have read: "avg = avg(n+1 : _______________________________________________ Help-octave mailing list [hidden email] https://www-old.cae.wisc.edu/mailman/listinfo/help-octave |
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
_______________________________________________ Help-octave mailing list [hidden email] https://www-old.cae.wisc.edu/mailman/listinfo/help-octave |
Powered by Nabble | Edit this page |