Polynomial matte matrix with some restrictions on the coefficients

Question

I have data that I must interpolate using a function, which should look like this:

f(x) = ax 4 + bx 2 + c

with a > 0 and b ≤ 0 . Unfortunately, MATLAB polyfit does not allow restrictions on polynomial coefficients. Does anyone know if there is a MATLAB function for this? Otherwise, how to implement it?

Thank you very much,

Elizabeth

Answer 1

You can use fminsearch , fminunc , defining your target function manually.

Alternatively, you can define your problem in a slightly different way:

f(x) = a 2 x 4 - b 2 x 2 + c

Now the new a and b can be optimized for constraints without , ensuring that the final a and b you are looking for are positive (negative or negative).).

Answer 2

Without limits, the problem can be written and solved as a simple linear system:

 % 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, [abc], found by % solving A*[abc]' = y in a least-squares sense abc = A\your_Y_data(:) 

These restrictions, of course, will be automatically satisfied iff that the limited model really underlies your data. For instance,

 % 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, [abc]) + 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)') 

However, if you place restrictions on data that is not exactly described by this limited model, things may go wrong:

 % 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, [abc]) + 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)') 

(this solution has a == 0 , which indicates the wrong model choice).

If exact equality a == 0 is a problem: it makes no difference if you set a == eps(0) . Quantitatively, this will not be noticeable for real data, but nonetheless it is non-zero.

In any case, I have a suspicion that your model is not selected well, and the limitations are a “correction” to make everything work, or your data must be objective or scaled before trying to adapt, or that some similar prerequisites apply ( I often saw people doing such things, so yes, I'm a little biased in this regard :).

So ... what are the real reasons for these restrictions?

Answer 3

If you have a toolbar for adjusting the curve, then fit allows you to set limits using the options 'upper' and 'lower'. You would like something like that.

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

Note the use of -inf to set a specific factor without restriction.

This will give a cfit object with appropriate coefficients. You can access them using, for example, M.p1 for the member x ^ 4. Alternatively, you can evaluate the function at any points that you want to use feval .

I think you can do a similar thing using lsqcurvefit in the optimization toolbar.