MATLAB >> using FMINSEARCH for weighted curve fitting

by ross » Sat, 09 Dec 2006 09:45:22 GMT

I asked a while ago about polynomial curve fitting using arbitrary weighting
functions for the fitted deviates.

John d'Errico came to my rescue with FMINSEARCH. Eg, using ERF as the
weighting function, we can search for the minimal weighted sums of squares as
follows:

% perform ordinary LS fit to get a starting value for
% the coeff vector for an n-th order poly
p = polyfit(x, n);

% define function to minimise
obj = @(c) sum( erf( (y - polyval(c,x)).^2 ) );

% perform simplex hull search for optimised coeff
% vector that gives local (n+1)D minimum in obj()
optim_p = fminsearch(obj, p);


In practice the spread of the weighting function needs to be matched to the
data. The simplest way to parameterise this is to scale the deviates before
weighting them. The definition of obj() is then:

obj = @(c) sum( erf( k * (y - polyval(c, x)).^2 ) );

Now, I don't know how the Nelder-Mead algorithm used by FMINSEARCH works, but I
gather the algorithm is meant to work generically, ie, the local minimum is
found by searching numerically in each dimension, without regard to what role
the variables play in the obj() function.

Thinking along these lines led me to combine the scaling factor k with the
polynomial coefficients into a composite vector to be optimised:
obj = @(c) sum( erf( c(1) * (y - polyval(c(2:end), x)).^2 ) );
optim_c = fminsearch(obj, [k p]);

If this is successful, then:
optim_k = optim_c(1);
optim_p = optim_c(2:end);

However, it ISN'T successful. The coeffs obtained are wildy wrong. Why
doesn't it work?


Thanks

Ross


MATLAB >> using FMINSEARCH for weighted curve fitting

by Dmitrey » Sat, 09 Dec 2006 15:41:54 GMT


Hi ross,
try OpenOpt - it contains ralg() - fminsearch analog; in many cases
it's much more better.
so, you can use
prob = ooAssign(@objFun,x0);
prob.useScaling = 1;
prob.ScaleFactor = (scale factor or, if unset, typicalX will be used,
or, if unset, x0);

r = ooRun(prob,ralg)

BTW it can use user-supplied (sub)gradient. In future I intend to add
linear & non-lin constraints, currently only lb-ub are available via
a wrapper.

currently SolvOpt has more powerfull implementation of ralg (with
preventing matrix B to be singular, and handling more constraints
than I), you can try that one (but it has very unconvinient
interface)

Latest OpenOpt ver 0.17 you can download from < http://www.box.net/public/6bsuq765t4> ;

However, if your funcs are smooth, other MATLAB funcs could be
better, for example lsqurvefit or nlinfit

WBR, D

MATLAB >> using FMINSEARCH for weighted curve fitting

by John D'Errico » Sat, 09 Dec 2006 18:50:37 GMT


I can tell you why it fails. In the
objective function expression

obj = @(c) sum( erf( c(1) * (y - polyval(c(2:end), x)).^2 ) );

drive c(1) to zero. What happens? Are
the coefficients of the polynomial even
relevant? (No.)

You can scale the problem this way, but
the scaling must be fixed in advance.

If you wish, why not use that first
polyfit call to get an initial set of
residuals? Then choose some value of k
that is reasonable.

A simple variation of this is just to
use a trimmed regression. Delete the
data points with the worst residuals
from the analysis. Drop the worst 25%
of the points from the analysis. Or,
if you are using a traditional weighted
regression, give them a zero weight.

John


--
The best material model of a cat is another, or preferably the same, cat.
A. Rosenblueth, Philosophy of Science, 1945

Those who can't laugh at themselves leave the job to others.
Anonymous

MATLAB >> using FMINSEARCH for weighted curve fitting

by ross » Sat, 09 Dec 2006 20:06:10 GMT


Thanks for your other words of advice, John.

But with respect to the explanation above, I don't understand.
In a few test runs, I found that c(1) is not returned as 0, but typically as a
number which is not in any way extreme, eg, 10-20.

I have been doing test runs with a setup similar to the one you mentioned, ie,
the data are produced by a simple exponential generator with added noise, the
fitting model is a 4th order poly, and the weighting function is ERF. I'm no
mathematician, but I think these specifications imply that the minimisation
surface is likely to be smooth and well behaved in the search region.
Therefore, if the scaling value c(1) has not run off to 0 or infinity, then why
shouldn't the search in the remaining dimensions (representing the poly coeffs)
settle smoothly on to a local minimum of obj() ? Furthermore, I am using good
starting estimates, so the local minimum should already be within sight, so to
speak, and the search algorithm shouldn't have to cross any pathological
terrain to reach the minimum.

cheers

ross

MATLAB >> using FMINSEARCH for weighted curve fitting

by John D'Errico » Sun, 10 Dec 2006 00:00:25 GMT


Actually, I think my original formulation
had the square outside the erf. It may be
less susceptible to numerical problems.




Large numbers for k are just as bad.

If you multiply your residuals by a
large enough number, erf may be
overwhelmed. Erf will return uniformly
the number 1.

erf(20)
ans =
1




It does not need to go to inf. 20 may
be entirely large enough to kill off
the response.



No. If k is large enough, changes in
the poly do not effectively change
the transformed residuals at all. Erf
is maxed out, so all of the residuals
are mapped to 1.



See above.

John


--
The best material model of a cat is another, or preferably the same, cat.
A. Rosenblueth, Philosophy of Science, 1945

Those who can't laugh at themselves leave the job to others.
Anonymous

Similar Threads

1. Curve Fitting tool: robust fit: view final weights

Hello,

I'm new to matlab, but managed to generate an M-code using the curve fitting 
tool box and even to modify this code.

What I couldn't find in the manual/online help is the syntax to display the final 
weights that the fit command uses in the robust mode. It would be nice to 
actually see or check how each data point is weighted in the regression.

Thanks for any help.
thorsten

Here are some code snipplets of the M-file I use:
fo_ = 
fitoptions('method','NonlinearLeastSquares','Robust','On','MaxFunEvals',15000);
ok_ = isfinite(x) & isfinite(y);

% Fit this model using new data
[cf_,gof] =  fit(x(ok_),y(ok_),ft_,fo_);


2. (curve fitting) which curve model is good for fitting this curve(image attached: 44KB) - MATLAB

3. Weight-bug in curve fitting toolbox?

Dear readers,
I have recently been looking through the options in Matlab R2008a for performing a weighted robust linear regression. It seems that the curve fitting toolbox is the only function that accepts additional weights. However, I might have come across some inconsistency in the usage of the weights that I would like to have confirmed.

In the help section on Robust Least Squares (Least-Squares Fitting, Curve Fitting Toolbox) it says the following: "Note that if you supply your own regression weight vector, the final weight is the product of the robust weight and the regression weight".

In the code of fit.m there is no such multiplication, in fact the weights from the robust fit are not eve given as outputs. What I did find though is that the initially supplied weights are squared on the first line in the "goodstruct"-function. Is this really correct or is this supposed to be the product of the robust and initial weights as stated in the help section? Unfortunately this calculation will affect all the subsequent calculations of the statistics so at the moment I don't really trust in using the cft for this purpose.

I have checked the code in Matlab 2009a as well and it is the same thing. I would really appreciate it if someone could tell me if this is correct or if this is a bug that has slipped through every revision. Thanks!

4. curve fitting with weights - MATLAB

5. Weight in curve fitting

Hi, does somebody knows how to put weight in just one point or two points in the data when iam trying to fit a curve ??? (using the editor, with commands, not the curve fitting tool ), thanks a lot!

6. Weights in curve fitting - MATLAB

7. curve fitting and f-values using fit()

Hi list,

Is there a way of getting an f-value or some form of significance
value on curves fitted with fit(). I've been looking but I don't seem
to find anything close to a standard f or p-value. The r-square is
available but this only tells me how close the values are to the
fitted curve not necessarily something about the strength of the fit.

Kind regards,
Koen

8. curve fit and add fitted curve together - MATLAB