MATLAB Answers

How to recreate a piecewise polynomial from polyfit to values in each interval?

68 views (last 30 days)
Will
Will on 15 Aug 2018
Commented: dpb on 16 Aug 2018
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'})

Accepted Answer

dpb
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).
  2 Comments
dpb
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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by