Re: leasqrdemo

[Michael Creel]

Przemek Klosowski wrote:
>    I am trying to fit function b*(x^a), but I don't know much about fitting
>    and I also have problems with derivatives.
>    I am using leasqrdemo.m, where I defined
>
>    p = [-1,10];
>    ....
>    function y = leasqrfunc(x,p)
>            y = p(2)*(x^p(1))
>    endfunction
>
>    but I don't know what the leasqrdfdp should look like.
>    function y = leasqrdfdp(x,f,p,dp,func)
>            y = [?,?]
>    endfunction
>
> You don't need to define the derivatives---if you don't they are
> calculated numerically. leasqr() is a wonderful facility, but it like
> a race car: temperamental and demanding.  Still, with little care, it
> is a pleasure to drive. Here's a simple example:
>
>     x=1:10; y=[12,5,3,2,1,1,1,1,1,1]; p = [-1,10];
>     function y = leasqrfunc(x,p)
>        y = p(2)*(x.^p(1))
>     endfunction
>
>     [f,pout]=leasqr(x',y',p,"leasqrfunc");
>
> results in a reasonable fit:    pout = [ -1.2841; 12.0054 ]
>
> Note that I changed your function slightly by using the 'scalar'
> exponentiation .^ rather than matrix exponentiation, just so that
> I can use the function leasqrfunc() on a vector of values, e.g.:
>
>   plot(x,y,x,leasqrfunc(x,pout),x,f)
>
> Here's my list of issues with leasqr():
>
>  - it seems to be unnecessarily picky about its input arguments: x,y 
>    apparently must be column vectors, as leasqr(x,y,p,"leasqrfunc")
>    returns an error.
>
>  - I can't seem to get it to run quietly: it produces output even
>    with a semicolon at the end of the command line.
>
>  - the help text is exhaustive but not beginner-friendly; leasqrdemo()
>    didn't help me either
>
> I think that leasqr() is an important Octave facility, and I would be
> glad to do some work on it. Is there someone else who likes it too, to
> discuss what needs done and how? I haven't used leasqr() that much, so
> maybe I am missing something.


leasqr implements a specific method of finding nonlinear least squares 
estimates. There are alternatives, for example, using BFGS to do the 
minimization. You just need to define a function that calculates the 
sum of squared errors rather than the fitted values. Code for the 
above example is
------------------------------------------------
x=1:10;
y=[12,5,3,2,1,1,1,1,1,1];

function y = leasqrfunc(x,p)
        y = p(2)*(x.^p(1));
endfunction


function obj_value = obj(p, y, x);
        fit = leasqrfunc(x,p);
        e = y - fit;
        obj_value = sumsq(e)/rows(e);
endfunction

# plain NLS fit using bfgsmin
p = zeros(2,1);
p = bfgsmin('obj', {p, y, x});
p
obj(p,y,x)
-------------------------------------------------

The output is
octave:1> example
p =

   -1.2842
   12.0054

ans = 0.57346
octave:2>

so the same result is obtained. The bfgsmin routine will let you use 
all zeros as start values. The reason leasqr was not running quietly 
is that there was a missing semicolon in the leasqrfunc().

Cheers, M.