Moving polynomial fit

[babelproofreader]

I would like to write a script/function for use on a time series that will least squares fit a polynomial of a given degree to a moving window across the time series in the same way as a moving average computes an average of a moving window, the output of the script/function being a vector of the extreme right hand edge/last values of the fitted polynomial over the moving window. I have tried using the Savitsky-Golay filters from the Signals package but they are not suitable for my purposes as the end conditions mean that the values of the fitted polynomial at the right hand edges of the moving window are different from the values calculated when the data are no longer at the right hand edge. What in-built Octave functions should I consider using to achieve the above i.e. polyfit or some other function(s)? Also, as loops can be quite slow, is there a more efficient way of coding the above rather than having to resort to a for loop?

[macy-3]

I had to do something like this for 2 minutes of streaming data once.  The
FOR loop took 20 minutes+ and the array process took 12 seconds!

It went something like this:

Multiply the original vector in such a way as to make a 'shifted' array.

Where each line of the array is your original vector, but shifted slightly.

Then, apply to the array in a single stomp and you should be there

I had to expand the memory octave used, else kept crashing my system when
the data files were up around 20MB

Robert

> I would like to write a script/function for use on a time series that will
> least squares fit a polynomial of a given degree to a moving window across
> the time series in the same way as a moving average computes an average of
> a
> moving window, the output of the script/function being a vector of the
> extreme right hand edge/last values of the fitted polynomial over the
> moving
> window. I have tried using the Savitsky-Golay filters from the Signals
> package but they are not suitable for my purposes as the end conditions
> mean
> that the values of the fitted polynomial at the right hand edges of the
> moving window are different from the values calculated when the data are
> no
> longer at the right hand edge. What in-built Octave functions should I
> consider using to achieve the above i.e. polyfit or some other
> function(s)?
> Also, as loops can be quite slow, is there a more efficient way of coding
> the above rather than having to resort to a for loop?
>
>
> --
> View this message in context:
> http://www.nabble.com/Moving-polynomial-fit-tp25309400p25309400.html
> Sent from the Octave - General mailing list archive at Nabble.com.

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Ben Abbott]

On Sep 5, 2009, at 11:09 AM, babelproofreader wrote:

> I would like to write a script/function for use on a time series  
> that will
> least squares fit a polynomial of a given degree to a moving window  
> across
> the time series in the same way as a moving average computes an  
> average of a
> moving window, the output of the script/function being a vector of the
> extreme right hand edge/last values of the fitted polynomial over  
> the moving
> window. I have tried using the Savitsky-Golay filters from the Signals
> package but they are not suitable for my purposes as the end  
> conditions mean
> that the values of the fitted polynomial at the right hand edges of  
> the
> moving window are different from the values calculated when the data  
> are no
> longer at the right hand edge. What in-built Octave functions should I
> consider using to achieve the above i.e. polyfit or some other  
> function(s)?
> Also, as loops can be quite slow, is there a more efficient way of  
> coding
> the above rather than having to resort to a for loop?

There might be a better way ... Assuming you'd like to center your  
window about the data point of interest ...

---------------- begin: mwpolysmooth.m ---------------
function ym = mwpolysmooth (x, y, order, window)
% USAGE: ym = mwpolysmooth (x, y, order, window)

   xc = num2cell(x);
   idx = @(xc) find (abs (x - xc) < 0.5*window);
   n = cellfun (idx, xc, 'uniformoutput', false);
   mpfit = @(n) polyfit (x(n), y(n), order);
   pm = cellfun (mpfit, n, 'uniformoutput', false);
   nc = num2cell (1:numel(x));
   mpval = @(n) polyval (pm{n}, x(n));
   ym = cell2mat (cellfun (mpval, nc, 'uniformoutput', false));

end

%!demo
%! x = 0:0.1:10;
%! y = rand (size(x)) + cos (x);
%! window = 5;
%! order = 2;
%! ym = mwpolysmooth (x, y, order, window);
%! clf
%! plot (x, y, 's', x, ym)
---------------- end: mwpolysmooth.m ---------------

Ben




_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[babelproofreader]

On Sep 5, 2009, at 11:09 AM, babelproofreader wrote:

> I would like to write a script/function for use on a time series  
> that will
> least squares fit a polynomial of a given degree to a moving window  
> across
> the time series in the same way as a moving average computes an  
> average of a
> moving window.....

There might be a better way ... Assuming you'd like to center your  
window about the data point of interest ...

---------------- begin: mwpolysmooth.m ---------------
function ym = mwpolysmooth (x, y, order, window)
% USAGE: ym = mwpolysmooth (x, y, order, window)

   xc = num2cell(x);
   idx = @(xc) find (abs (x - xc) < 0.5*window);
   n = cellfun (idx, xc, 'uniformoutput', false);
   mpfit = @(n) polyfit (x(n), y(n), order);
   pm = cellfun (mpfit, n, 'uniformoutput', false);
   nc = num2cell (1:numel(x));
   mpval = @(n) polyval (pm{n}, x(n));
   ym = cell2mat (cellfun (mpval, nc, 'uniformoutput', false));

end
---------------- end: mwpolysmooth.m ---------------

Ben,

Thanks for your function - I've learned some useful coding from it - but I'm not sure if it really address my problem. I wish to fit real time streaming data (market prices), but as I understand it your code "peaks into the future" to smooth over the window i.e. the fitted value for the most recent streaming data "n" is different from the value that will be calculated when this data value is (n - some_value) in the past, hence invalidating any historical testing over past prices. The simple for loop below is an example of what I want to do, but it is slow.

x=(1:5)';
polysmooth=price;
for i=5:length(price)
y=price(i-4:i,1);
a=polyfit(x,y,2);
polysmooth(i,1)=polyval(a,5);
endfor

[macy-3]

NOW I REMEMBER!

I had to rewrite the function polyval into a custom type function so that
I could use c variable as an array, not just a vector.

Then, just write the whole for-loop as an array sequence.  As I said I
went from 20 minutes down to 12 seconds for mine.

Robert

> On Sep 5, 2009, at 11:09 AM, babelproofreader wrote:
>
>> I would like to write a script/function for use on a time series
>> that will
>> least squares fit a polynomial of a given degree to a moving window
>> across
>> the time series in the same way as a moving average computes an
>> average of a
>> moving window.....
>
> There might be a better way ... Assuming you'd like to center your
> window about the data point of interest ...
>
> ---------------- begin: mwpolysmooth.m ---------------
> function ym = mwpolysmooth (x, y, order, window)
> % USAGE: ym = mwpolysmooth (x, y, order, window)
>
>    xc = num2cell(x);
>    idx = @(xc) find (abs (x - xc) < 0.5*window);
>    n = cellfun (idx, xc, 'uniformoutput', false);
>    mpfit = @(n) polyfit (x(n), y(n), order);
>    pm = cellfun (mpfit, n, 'uniformoutput', false);
>    nc = num2cell (1:numel(x));
>    mpval = @(n) polyval (pm{n}, x(n));
>    ym = cell2mat (cellfun (mpval, nc, 'uniformoutput', false));
>
> end
> ---------------- end: mwpolysmooth.m ---------------
>
> Ben,
>
> Thanks for your function - I've learned some useful coding from it - but
> I'm
> not sure if it really address my problem. I wish to fit real time
> streaming
> data (market prices), but as I understand it your code "peaks into the
> future" to smooth over the window i.e. the fitted value for the most
> recent
> streaming data "n" is different from the value that will be calculated
> when
> this data value is (n - some_value) in the past, hence invalidating any
> historical testing over past prices. The simple for loop below is an
> example
> of what I want to do, but it is slow.
>
> x=(1:length(price))';
> polysmooth=price;
> for i=5:length(price)
> y=price(i-4:i,1);
> a=polyfit(x,y,2);
> polysmooth(i,1)=polyval(a,5);
> endfor
>
>

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Ben Abbott]

On Sep 7, 2009, at 10:27 AM, babelproofreader wrote:

>
>
> On Sep 5, 2009, at 11:09 AM, babelproofreader wrote:
>
>> I would like to write a script/function for use on a time series
>> that will
>> least squares fit a polynomial of a given degree to a moving window
>> across
>> the time series in the same way as a moving average computes an
>> average of a
>> moving window.....
>
> There might be a better way ... Assuming you'd like to center your
> window about the data point of interest ...
>
> ---------------- begin: mwpolysmooth.m ---------------
> function ym = mwpolysmooth (x, y, order, window)
> % USAGE: ym = mwpolysmooth (x, y, order, window)
>
>   xc = num2cell(x);
>   idx = @(xc) find (abs (x - xc) < 0.5*window);
>   n = cellfun (idx, xc, 'uniformoutput', false);
>   mpfit = @(n) polyfit (x(n), y(n), order);
>   pm = cellfun (mpfit, n, 'uniformoutput', false);
>   nc = num2cell (1:numel(x));
>   mpval = @(n) polyval (pm{n}, x(n));
>   ym = cell2mat (cellfun (mpval, nc, 'uniformoutput', false));
>
> end
> ---------------- end: mwpolysmooth.m ---------------
>
> Ben,
>
> Thanks for your function - I've learned some useful coding from it -  
> but I'm
> not sure if it really address my problem. I wish to fit real time  
> streaming
> data (market prices), but as I understand it your code "peaks into the
> future" to smooth over the window i.e. the fitted value for the most  
> recent
> streaming data "n" is different from the value that will be  
> calculated when
> this data value is (n - some_value) in the past, hence invalidating  
> any
> historical testing over past prices.

You can modify the window to represent only past values by ...

   idx = @(xc) find ((xc - x) < window & (xc - x) > 0);

Ben



_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave