Dopasowywanie wielomianu do maksimum funkcji

Question

Chcę dopasować wielomian do hałaśliwych danych, tak aby przybliżony wielomian był zawsze> = oryginalne dane. Na przykład:

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 

Czy istnieje lepszy sposób/metoda, która również pracuje z wyższymi wielomianów porządku?

Sposób easies byłoby po prostu użyć polyfit a następnie obliczyć max(y-yi) i dodać to jako przesunięcie, ale nie byłoby to optymalne ...

EDIT: Chcę używać GNU Octave, ale dodał „Matlab” jako tag, ponieważ język i funkcje są podobne.

EDIT: na podstawie odpowiedzi TVO i danych rzeczywistych:

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 

ale jak widać, zoptymalizowany i oceniane poli posiada punkty < punktowych orginal co jest niedozwolone.

Answer 1

Twoja metoda wydaje się być w porządku i nie widzę powodu, dla którego nie można jej użyć do wielomianów wyższego rzędu. Proszę wyjaśnić, dlaczego uważasz, że nie można go użyć.

Używasz narzędzia "sqp" programu Octave. Dokumentacja jest tutaj: http://www.gnu.org/software/octave/doc/v4.0.1/Nonlinear-Programming.html

Możesz uniknąć pomnożenia błędu przez dowolną liczbę (1000 w twoim przykładzie), gdy jest ujemna. To może się nie udać w przypadku różnych zestawów danych, zwłaszcza jeśli są one większe, tj. Więcej punktów danych.

Możesz spróbować użyć nieliniowych opcji ograniczenia nierówności, które oferuje "sqp" Octave'a, tj. H (x)> = 0 (patrz dokumentacja).

Jako funkcja celu phi można użyć błędu kwadratowej normy, jak w przykładzie, i dodać ograniczenia formularza h (x)> = 0 dla każdego punktu danych. Zauważ, że "x" byłoby współczynnikami wielomianowymi, które chcesz dopasować, a h (x) jest wielomianem ocenianym w określonym punkcie danych.

Na przykład:

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); 

Wskazówka że funkcja przymusu „h” zapewnia, że ​​wielomian będzie leżeć powyżej (x1, y1) i funkcja celu „phi” będzie starał się utrzymać go jak najbliżej możliwy. Możesz rozszerzyć "h" o jedno ograniczenie dla każdego punktu danych w swoim zestawie (patrz dokumentacja).