Re: proposal for new m-file function

[Olaf Till]

On Thu, Mar 22, 2012 at 10:06:52AM +0100, Olaf Till wrote:
> On Wed, Mar 21, 2012 at 08:33:01PM -0400, Ben Abbott wrote:
> > A user's question about fixed points piecewise-linear fitting led to the 
> > development of a function looks to be a good fit for Octave's core.
> > 
> >     
> > https://mailman.cae.wisc.edu/pipermail/help-octave/2012-March/050900.html
> > 
> > The idea was to fit a piece-wise polynomial to a set of data. After some 
> > discussion between myself and Martin Helm, the attached ppfit.m was 
> > produced. The name was chosen to match the existing ppval().
> > 
> > The function provides a least-squares fit of a 1D interpolation with 
> > specified break positions to a set of data.
> > 
> > Demos and tests are included.
> > 
> > Any concern about adding this to Octave's core ?
> > 
> > Ben
> 
> The loop could be vectorized:
> 
> --- original_ppfit.m    2012-03-22 09:18:22.000000000 +0100
> +++ ppfit.m     2012-03-22 09:32:04.000000000 +0100
> @@ -113,12 +113,8 @@
>             "Input data value pairs must have the same size.")
>    endif
>    ## Use linear least squares for the initial guess
> -  a = zeros (numel (x), numel (xi));
> -  for ii = 1:numel (xi)
> -    yi = zeros (numel (xi), 1);
> -    yi(ii) = 1;
> -    a(:, ii) = interp1 (xi, yi, x, method, "extrap");
> -  endfor
> +  yi = eye (numel (xi));
> +  a = interp1 (xi, yi, x(:), method, "extrap");
>    ## y = (q*r) * yi + errors
>    ws = warning ("off", "all"); 
>    [q, r, p] = qr (a .* w, 0);
> 
> but for a bug in interp1 which transposes the returned matrix if both
> "spline" and "extrap" are given.
> 
> Olaf

Attached is the bugfix for interp1.m.

Olaf

-- 
public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net

interp1.changeset
Description: Text document

signature.asc
Description: Digital signature