Re: NMinimize fails to find a minimum value
- To: mathgroup at smc.vnet.net
- Subject: [mg120859] Re: NMinimize fails to find a minimum value
- From: Ray Koopman <koopman at sfu.ca>
- Date: Sat, 13 Aug 2011 06:48:42 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- References: <j1im1r$cmu$1@smc.vnet.net>
On Aug 5, 11:13 pm, Hovhannes Babayan <bhovhan... at gmail.com> wrote:
> Hello,
>
> I am solving the following task. I have set of points (list of pair of
> X, Y coordinates) and list of weights for the each point. If for some
> point weight=1, that point is very important of weight=0, that point
> is not important. Weight ranges from 0 to 1.
> The task is to construct a curve, which will fit given list of points
> taking weights into account. The analytic form of function is known,
> but I have to find parameters z,u,a,b, which will give the best
> fitting curve.
> To do that, I am using NMinimize to minimize the following expression:
> Sum[ ((f[x[[i]],z,u,a,b]-y[[i]])^2)*weights[[i]], {i, n} ]
> where "n" is a number of points, "x" and "y" are point coordinates,
> and "f(r, z,u,a,b)" is my function, where z,u,a,b are parameters which
> should be found and "r" (x) is the only argument.
>
> [...]
>
> 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 *)
>
> [...]
I think you're making things more difficult than they need to be.
1. NMinimize doesn't need to know about z. Since it's just a
multiplier on the approximation and you're doing least-squares
fitting, there is a closed-form solution for z given a, b, and u.
2. The parameters z,u,a,b in f may make good substantive sense,
but d * r^a * (r+c)^b is algebraically equivalent to f and
is easier for NMinimize (and similar routines) to work with.
Again, there is a closed-form solution for d given a, b, and c.
3. Think in terms of vectors instead of scalars.
Use Dot to sum products.
Here's the code I used and the results I got.
The minimand is the square root of the weighted mean square error.
'e' is the effective number of points. An argument can be made
that multiplying the minimand by Sqrt[e/(e-k)], where k is
the number of parameters estimated -- in this case, 4 -- would
give a more honest appraisal of the size of the errors. That
adjustment would not change the parameter estimates, though.
The plot makes it clear that the last point is basically ignored,
which the very low weight assigned to it suggests might happen.
{x,y} = Transpose@points; w = weights/Tr@weights
e = 1/w.w
{.0472525, .0750459, .108645, .115987, .127318, .187287,
.149024, .0963404, .0693524, .0202083, .00353972}
8.25573
g[r_, a_,b_,c_] := r^a (r+c)^b;
Clear[a,b,c]; wy = w*y;
NMinimize[t = g[x,a,b,c]; Sqrt[w.(t*wy.t/w.t^2 - y)^2], {a,b,c}]
{a,b,c} = {a,b,c}/.%[[2]]; t = g[x,a,b,c]; d = wy.t/w.t^2
{.00338027, {a -> -.435854, b -> -2.44512, c -> 30.882}}
1880.97
Plot[g[r,a,b,c]*d,{r,4,56}, Prolog->MapThread[
{AbsolutePointSize@#1,Point@#2}&,{40*Sqrt@w,points}]]