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

curve fittingMATLABpiecewise-polynomialpolynomial fittingspline

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'})

Best Answer

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).
Related Question