Re: non-linear fitting with conditions

[Olaf Till]

On Mon, Mar 22, 2010 at 01:48:49AM -0400, Carnë Draug wrote:
> Hi
> 
> the function leasqr (part of the optim package) is used to fit non-linear
> data. In my case, some values for the parameters (there's 3 of them) make no
> physical sense so I made a condition inside the function so it returns
> before evaluating such values.
> 
> if (rCon < 0 || sigma < 0 || theta < 0 || theta  > 1)
>    return
> endif
> 
> However, when leasqr attempts to fit such values, it returns an error (at
> least that's my understanding of the origin of the error, when I looked into
> the source) which makes sense because he tries to subtract two values, one
> with zero size. Here's the error
> 
> error: operator -: nonconformant arguments (op1 is 1x26, op2 is 0x0)
> error: called from:
> error:   /usr/share/octave/packages/3.2/optim-1.0.10/leasqr.m at line 320,
> column 3
> 
> If I got it right, op1 = y (argument of leasqrand) and op2 = f (return of
> feval with function F).
> 
> I guess it makes sense for the code to check this but in that case how can I
> check for this and make sure that leasqr doesn't give me certain values? Or
> should this be considered a bug?
> I though about making a while loop with a try statement inside to catch any
> errors and change the initial guess of the parameters values until leasqr
> doesn't try one that makes no sense on its iterations. But I fear that
> sometimes this will create an infinite loop.
> 
> By the way, here's the function that I am using
> 
> function y = func_radial_profile(pf_distances, parameters)
> 
> rCon  = parameters(1);
> sigma = parameters(2);
> theta = parameters(3);
> 
> if (rCon < 0 || sigma < 0 || theta < 0 || theta  > 1)
>     return
> endif
> 
> for i=rows(pf_distance)
>     y(i) = 1-(1-theta)*exp(-(pf_distances(i)-rCon)^2/(2*sigma^2))
> endfor

A usual way to cope with such problems is setting bounds for
parameters. In leasqr, this is done in your case with:

options.bounds = [0, Inf; \ # rCon
                  0, Inf; \ # sigma
                  0, 1]; # theta

and calling leasqr with "options" as 10th argument:

[..., ...] = leasqr (x, y, pin, "func_radial_profile", .0001, \
                     20, ones (3, 1), .001 * ones (3, 1), "dfdp", \
                     options);

However, you should use leasqr of package optim-1.0.12 for that, since
leasqr of optim-1.0.10 will in general not converge if bounds are set
(one might not even notice this, since the result might look good,
although it isn't).

Olaf