Re: Moving polynomial fit

[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