Set polynomial with function maxima

Question

I want to fit the polynomial to noisy data, so that the approximated polynomial is always> = source data. For example:

x = linspace (-2, 6); y = (x-2).^2 + 1 + 2 * randn (size (x)); function ret = delta (P, x, y) yP = polyval (P, x); d = yP - y; d (d < 0) *= 1000; ret = sumsq (d); endfunction P0 = polyfit (x, y, 2); f = @(P) delta (P, x, y); [P, FVAL] = sqp (P0, f) xi = linspace (min(x), max(x), 100); yi = polyval (P, xi); plot(x, y, xi, yi); grid on 

Is there a better way / method that also works with higher order polynomials?

An easy way would be to simply use polyfit and then calculate max(y-yi) and add this as an offset, but that would not be optimal ...

EDIT: I want to use GNU OCtave, but added "matlab" as a tag because the language and functions are similar.

EDIT: based on tvo response and real data:

 x = [10 20 30 40 50 60 80 100]; y = [0.2372, 0.1312, 0.0936, 0.0805, 0.0614, 0.0512, 0.0554, 0.1407]; function ret = delta (P, x, y) ret = sumsq (polyval (P, x) - y); endfunction f = @(P) delta (P, x, y); h = @(P) polyval(P, x) - y; P0 = polyfit (x, y, 3); [P] = sqp (P0, f, [], h) xi = linspace (min(x), max(x)); yi = polyval (P0, xi); yio = polyval (P, xi); plot(x, y, xi, yi, ";initial guess;", xi, yio, ";optimized;"); grid on 

but, as you can see, the optimized and rated poly has points <the original point, which is not allowed.

Answer 1

your method seems wonderful, and I see no reason why it cannot be used for higher order polynomials. Please explain why, if you think it should not be used.

You are using the Octave 'sqp' solver. The documentation is here: http://www.gnu.org/software/octave/doc/v4.0.1/Nonlinear-Programming.html

You might want to avoid multiplying the error by an arbitrary number (1000 in your example) when it is negative. This may fail for different data sets, especially if they are larger, that is, more data points.

You can try to use the nonlinear inequality constraint options that Octave 'sqp' suggests, i.e. h (x)> = 0 (see doc).

You can use the square norm error as the target function phi, as in your example, and add constraints of the form h (x)> = 0 for each data point. Note that “x” is the polynomial coefficients you want to fit, and h (x) is the polynomial evaluated at a particular data point.

For example:

 phi = @(P) delta_mod (P, x, y); % mod: Don't increase the importance of negative residuals!! h = @(P) polyval(P, x1) - y1; Psol = sqp(P0, phi, [], h); 

Note that the constraint function "h" ensures that the polynomial is higher than (x1, y1), and the objective function "phi" will try to keep it as close as possible. You can extend 'h' to contain one constraint for each data point in your set (see Doc).