Polynomial MATLAB forma con alcuni vincoli sui coefficienti

Question

Ho dati che devo interpolare con una funzione che deve essere del tipo seguente:

f(x) = ax4 + bx2 + c

con a > 0 e b ≤ 0. Sfortunatamente, polyfit di MATLAB non consente vincoli sui coefficienti del polinomio. Qualcuno sa se c'è una funzione MATLAB per fare questo? Altrimenti, come posso implementarlo?

Grazie mille in anticipo,

Elisabetta

Answer 1

Puoi provare a utilizzare fminsearch, fminunc definire la funzione obiettivo manualmente.

In alternativa, è possibile definire il problema leggermente diverso:

f(x) = a2x4 - b2x2 + c

Ora, il nuovo a e b possono essere ottimizzati per senza vincoli, garantendo nel contempo che la finale a e b che stai cercando sono positivi (resp. Negativo).

Answer 2

senza vincoli, il problema può essere scritto e risolto come un sistema lineare semplice:

% Your design matrix ([4 2 0] are the powers of the polynomial) 
A = bsxfun(@power, your_X_data(:), [4 2 0]); 

% Best estimate for the coefficients, [a b c], found by 
% solving A*[a b c]' = y in a least-squares sense 
abc = A\your_Y_data(:) 

Quei vincoli saranno ovviamente automaticamente soddisfatta se e solo se che costretto modello alla base infatti i vostri dati. Ad esempio,

% some example factors 
a = +23.9; 
b = -15.75; 
c = 4; 

% Your model 
f = @(x, F) F(1)*x.^4 + F(2)*x.^2 + F(3); 

% generate some noisy XY data 
x = -1:0.01:1; 
y = f(x, [a b c]) + randn(size(x)); 

% Best unconstrained estimate a, b and c from the data 
A = bsxfun(@power, x(:), [4 2 0]); 
abc = A\y(:); 

% Plot results 
plot(x,y, 'b'), hold on 
plot(x, f(x, abc), 'r') 
xlabel('x (nodes)'), ylabel('y (data)') 

Tuttavia, se si impongono vincoli su dati che sono non accuratamente descritto da quel modello vincolato, le cose potrebbero andare male:

% Note: same data, but flipped signs 
a = -23.9; 
b = +15.75; 
c = 4; 

f = @(x, F) F(1)*x.^4 + F(2)*x.^2 + F(3); 

% generate some noisy XY data 
x = -1:0.01:1; 
y = f(x, [a b c]) + randn(size(x)); 

% Estimate a, b and c from the data, Forcing a>0 and b<0 
abc = fmincon(@(Y) sum((f(x,Y)-y).^2), [0 0 0], [-1 0 0; 0 +1 0; 0 0 0], zeros(3,1)); 

% Plot results 
plot(x,y, 'b'), hold on 
plot(x, f(x, abc), 'r') 
xlabel('x (nodes)'), ylabel('y (data)') 

(questa soluzione ha a == 0, indicativo di un scelta del modello errato).

Se l'uguaglianza esatta di a == 0 è un problema: non c'è ovviamente alcuna differenza se si imposta a == eps(0). Numericamente, questo non sarà evidente per i dati del mondo reale, ma è comunque diverso da zero.

In ogni caso, ho il sospetto che il tuo modello non sia ben scelto e che i vincoli siano una "correzione" per far funzionare tutto, o che i tuoi dati dovrebbero essere effettivamente imparziale/riscalati prima di provare a fare qualsiasi adattamento, o che alcuni si applicano precondizioni simili (ho visto spesso persone fare questo genere di cose, quindi sì, sono un po 'prevenuto in questo senso :).

Quindi ... quali sono le vere ragioni dietro quei vincoli?

Answer 3

Se si dispone della casella degli strumenti di adattamento della curva, quindi fit consente di impostare i vincoli utilizzando le opzioni "superiore" e "inferiore". Vorresti qualcosa di simile.

M=fit(x, f, 'poly4', 'upper', [-inf, 0, -inf, 0, -inf], 'lower', [0, 0, 0, 0, -inf]); 

Nota usa -inf per impostare un coefficiente particolare da non vincolare.

Questo darà un oggetto cfit con i coefficienti pertinenti. È possibile accedere a questi utilizzando per esempio M.p1 per il termine x^4. In alternativa è possibile valutare la funzione in qualsiasi punto si desideri utilizzando feval.

Penso che si possa fare una cosa simile usando lsqcurvefit anche nella toolbox di ottimizzazione.