Re: Strange behaviour of NMinimize when doing curve fitting
- To: mathgroup at smc.vnet.net
- Subject: [mg120727] Re: Strange behaviour of NMinimize when doing curve fitting
- From: "Oleksandr Rasputinov" <oleksandr_rasputinov at hmamail.com>
- Date: Sun, 7 Aug 2011 06:16:59 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- References: <j1im4n$cni$1@smc.vnet.net>
On Sat, 06 Aug 2011 07:14:47 +0100, Hovhannes Babayan
<bhovhannes at gmail.com> wrote:
> Hello,
>
> I want to fit a set of points with some curve, which form I know, but
> coefficients have to be determined to make the best possible fit. Also,
> for the each point a weight (a number ranging from 0 to 1) is defined.
> Point with greater weight is more important and point with lower weight
> is not so important. Point with weight=0 doesn't matter at all, point
> with weight=1 is extremely important.
> The code below uses NMinimize and draw the best fit graphic.
> The function f[r_, z_, u_, a_, b_] is a fitting function, where z,u,a,b
> are parameters.
> The question is following. Notice that I've used a/110 and b/100,
> instead of just a and b. When I use just a and b, produced curve doesn't
> fit at all my points. When I use a/50 and b/50 result is better, with
> a/70 and b/70 even better, and so on. The best seems to be a/110 and
> b/100 combination. Why it is so important for NMinimize??
>> From theoretical point of view it absolutely doesn't matter what is
>> used: a/100 or a/1000 or just a. Parameter "a" should be found in every
>> case.
> Please, if you have met such behavior, let me know, because currently I
> am doing NMinimize of NMinimize just to find good values to divide a and
> b, and it has no sense :)
> In order to try the code below, just copy/paste it into Mathematica, it
> should print best fit result (first element is error, second element is
> array with values of parameters z,u,a,b), should plot points and fitting
> curve. Use a/1 and b/1 to have the worst result and a/110 and b/100 to
> have good result.
> Thanks forehand for any suggestions, remarks and/or ideas.
>
>
>
> points = {
> { 4.68, 0.156},
> { 7.63, 0.100},
> {10.56, 0.074},
> {13.48, 0.061},
> {16.55, 0.048},
> {20.33, 0.029},
> {25.39, 0.023},
> {30.28, 0.018},
> {36.05, 0.016},
> {42.60, 0.016},
> {55.49, 0.028}
> } ;
> weights = {0.2523, 0.4007, 0.5801, 0.6193, 0.6798, 1, 0.7957, 0.5144,
> 0.37030, 0.1079, 0.0189};
>
> f[r_, z_, u_, a_, b_] :=
> z*(r^(-a/110)*(1 + r/u)^(b/100)); (* Here z,u,a,b are parameters *)
>
> x = points[[All, 1]];
> y = points[[All, 2]];
> n = Length[points];
> bestFitVals = NMinimize[
> Sum[
> ((f[x[[i]], z, u, a, b] - y[[i]])^2)*weights[[i]],
> {i, n}
> ],
> {z, u, a, b}
> ]
> Show[
> Plot[
> Tooltip[f[r, z, u, a, b] /. bestFitVals[[2]]],
> {r, Max[xdata], Max[ydata]}
> ],
> ListPlot[points],
> PlotRange -> {{0, Automatic}, {0, 0.2}},
> AxesLabel -> {"r", "f[r]"}
> ]
>
Such behaviour is normal with numerical global optimization; the choice of
objective function is usually very important, especially in cases such as
this one where the values of the parameters are strongly correlated with
one another. However, at least part of the problem also comes from your
function being complex for certain values of the parameters. If you
instead minimize the residuals you get by taking only the real part of
your function and add the constraint that its imaginary part should be
equal to zero, you will obtain much better results.
Also, some numerical global optimization methods perform better than
others depending on the application. As it happens, the differential
evolution method was invented specifically for determining parameters of
strongly nonlinear functions subject to constraints on the function
values, and as such is very appropriate for your problem.
Putting it all together (and using your points and weights), we get:
f[r_, z_, u_, a_, b_] := z*(r^(-a)*(1 + r/u)^b);
x = points[[All, 1]];
y = points[[All, 2]];
n = Length[points];
bestFitVals = NMinimize[
{
Sum[((Re@f[x[[i]], z, u, a, b] - y[[i]])^2)*weights[[i]], {i, n}],
Sum[Im@f[x[[i]], z, u, a, b]^2, {i, n}] == 0
}, {z, u, a, b}, Method -> "DifferentialEvolution"
]
Show[
Plot[
Tooltip[f[r, z, u, a, b] /. bestFitVals[[2]]],
{r, Max[x], Max[y]}
],
ListPlot[points],
PlotRange -> {{0, Automatic}, {0, 0.2}},
AxesLabel -> {"r", "f[r]"}
]