68 views (last 30 days)

Show older comments

I need to recreate the piecewise polynomial coefficients from values of the evaluated piecewise polynomial. I want to create an equivalent pp struct.

The pp has been stored as values sampled to enable a unique fit of each polynomial segment (i.e. values at the break times, and 4 internal points for each 6th order polynomial segment: 6 values, 6 polynomial coefficients to solve).

I thought I could use polyfit() on each break interval of samples to get each segment's polynomial coefficients, then use mkpp(). Scaled and centred, or not, the polyfit() coefficients in pp form do not reproduce the original samples, but evaluating polyval(polyfit()) for each segment does. I'm not clear why there's a difference.

What have I missed? Can it be done without the Curve Fitting Tool box?

An example for a single 6th order polynomial, split into 3 segments shows the discrepancy:

% Define polynomial coefs, x and y

p = [1, 2, 3, 4, 5, 6];

x = -5:10;

y = polyval(p, x);

% Define a pp form with knots at -5, 0, 5, 10

d = 1;

knotInterval = 5;

breaks = -5:knotInterval:10;

pieces = length(breaks)-1;

coefs = nan(pieces*d, length(p));

coefsScaled = nan(pieces*d, length(p));

mu = nan(pieces*d, 2);

yPolyval = nan(pieces, knotInterval+1);

% For the 6 (x,y) values in each pp break interval, fit an order 6

% polynomial, centered and scaled polynomial, and evaluate the standard

% fitted polynomial.

for nPiece = 1:pieces

idxPiece = x >= breaks(nPiece) & x <= breaks(nPiece+1);

coefs(nPiece, :) = polyfit(x(idxPiece), y(idxPiece), length(p)-1);

[coefsScaled(nPiece, :), ~, mu(nPiece, :)] = polyfit(x(idxPiece), y(idxPiece), length(p)-1);

yPolyval(nPiece,:) = polyval(coefs(nPiece, :), x(idxPiece));

end

% Build pp for the coefs and scaled coefs

pp = mkpp(breaks, coefs, d);

ppScaled = mkpp(breaks, coefsScaled, d);

% Plot the original data, and the fitted pp and polynomials

xPolyval = [-5:0; 0:5; 5:10];

plot(x, y, 'x-', x, ppval(pp, x), 'o--', x, ppval(ppScaled, x), '^--', ...

xPolyval(1, :), yPolyval(1, :), 'sq-.', xPolyval(2, :), yPolyval(2, :), 'sq-.', ...

xPolyval(3, :), yPolyval(3, :), 'sq-.')

ylim(1e5*[-1, 3])

legend({'True data', 'pp', 'ppScaled', ...

'Recovered polyval 1', 'Recovered polyval 2', 'Recovered polyval 3'})

dpb
on 16 Aug 2018

Edited: dpb
on 16 Aug 2018

OK, I did have time to go read ppval; indeed, it does presume coefficients are evaluated from each section as the origin, not in absolute coordinates--inside the function is

...

[b,c,l,o]=unmkpp(pp); % unpacks pp-->[bkpt,coefs,pieces,order]

[~,index] = histc(xs,[-inf,b(2:l),inf]); % get which section each input xq goes with

...

% now go to local coordinates ... % subtracts beginning breakpoint x from input

xs = xs-b(index); % vector of x so it thinks each is at origin

...

and goes on to evaluate the functional from there.

To make it work as you've defined the coefficients and coordinates, you have to evaluate over a modified xq vector that counteracts this internal shifting...

[~,ix]=histc(x,[-inf pp.breaks(2:pp.pieces) inf]); % get the equivalent index vector

xq=x+pp.breaks(ix); % "real" x plus breakpoint by section

with that then,

plot(x,y,'*-')

hold on

plot(x,ppval(pp,xq),'r-')

returns the expected result

Alternatively, compute the coefficients for each region based on the breakpoint for the region so each section origin is the breakpoint for the section.

Have to read the details of the documentation very carefully; the examples don't show that detail at all clearly. The description for the coefs though does spell it out; of course I didn't read that until just now after all the above! :)

coefs -- Poynomial coefficients

Polynomial coefficients, specified as an L-by-k matrix with the ith row coefs(i,:) containing

the local coefficients of an order k polynomial on the ith interval, [breaks(i), breaks(i+1)].

In other words, the polynomial is

coefs(i,1)*(X-breaks(i))^(k-1) + coefs(i,2)*(X-breaks(i))^(k-2) + ... + coefs(i,k-1)*(X-breaks(i)) + coefs(i,k).

dpb
on 16 Aug 2018

Yes, agreed. For higher-order polynomials conditioning is key.

It would take a modified ppval to use the centered form of coefficients in polyfit/val but that wouldn't be too difficult to do; it is actually very straight-ahead code once one looks at it; it clearly hasn't been infected by the later penchant to use all the fancy constructs and internals. :)

It might well be worth looking at it and "rolling your own" version if you're doing this regularly and with large x values in order to avoid numerical issues.

You could still use the structure form; there would be no issues I see with adding another field to the pp structure as flag; there's an internal logic branch now to treat returns from interp1 vis a vis mkpp and friends as far as the output orientation of results if the input xq values are array rather than vector; adding a custom field for 'centered' and the necessary statistics wouldn't interfere with that. You'd probably want to incorporate into a custom mkpp as well for completeness.

While old functions, seems like a worthwhile enhancement request on its own in general, not just for this type of use.

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!