Beefy Boxes and Bandwidth Generously Provided by pair Networks
Clear questions and runnable code
get the best and fastest answer
 
PerlMonks  

Re: perl minpack bindings or other function fitting code?

by swl (Parson)
on Nov 06, 2017 at 02:49 UTC ( [id://1202794]=note: print w/replies, xml ) Need Help??


in reply to perl minpack bindings or other function fitting code?

Given the lack of responses, I'll assume there are no bindings, or at least none that are obvious. Some further searching has turned up a few related modules, but none that implement minpack type functionality (see below).

In terms of the specific problem I described, the fitting issues were partly a PEBKAC, and partly the general approach to derivatives used in the various fit routines.

The PDL::Fit::LM code in my example needed a bounds constraint to ensure the partial derivative with respect to $c avoided log(0). Setting it to $c = max (0.000000001, $c) before any calculations resulted in convergence on the first attempt. I'd already tried bounds in an earlier iteration to avoid negative values (line was commented out in the code) but used 0 instead of a small positive value, hence the partial derivative with respect to $c could be set to infinity because that is how PDL handles log(0), as opposed to croaking as plain perl does.

PDL::Fit::Levmar also supports bounds, but I did not get it to work after a cursory attempt. The fact that it does not install cleanly via cpanm means that I won't look at it further for my needs in any case. search.cpan.org also cannot locate it, but it is at http://search.cpan.org/~jlapeyre/PDL-Fit-Levmar-0.0096/ if people want to see the details.

After continuing to look into the general MINPACK process, that library uses finite difference methods to calculate the derivatives instead of user specified functions. It should be possible to adapt the PDL routine to do the same, but that's probably something to take to the PDL list if I don't work it out myself.

And to provide a bit more detail in answer to my initial query, there are levenberg-marquardt fitting routines in the GNU Standard Library, but they are currently not included in the perl bindings. I don't know if there are plans to do so, but assume that if it were easy to do then it would likely have been done already.

The CEPHES library also includes the lmdif code from MINPACK, but the Math::CEPHES library does not include that as part of the subset it provides.

The final point perhaps worth noting is that the bounded PDL code takes 68 iterations to converge given my original starting parameters, whereas the minpack code via R takes 16. I'm assuming this is primarily due to differences in the settings used, for example the step size, but the minpack code does use other modifications. Setting the starting values to [-1,1,1,2.27395324135227] causes PDL::Fit::LM to use 22 iterations and the R version to use 21.

Replies are listed 'Best First'.
Re^2: perl minpack bindings or other function fitting code?
by swl (Parson) on Mar 25, 2018 at 22:42 UTC

    This is an update for people who find this thread in the future.

    PDL-Fit-Levmar has been updated and is now indexed on metacpan.

    It also works best of the modules I tried, converging on a solution without being sensitive to starting conditions.

    My previous issues were due to my own syntax errors. The main thing is that the code examples use the form x = f(t), not y = f(x) which I was assuming when I was testing, and hence my fits were using reversed axes. Such are the consequences of not reading the manual closely enough...

Re^2: perl minpack bindings or other function fitting code?
by etj (Deacon) on May 07, 2022 at 23:04 UTC
    We're always open to more/better functions in the GSL bindings for PDL! If you'd like to help out (I promise it's pretty easy, the existing bindings offer a great pattern to follow), please join one of the mailing lists - see PDL for links.

Log In?
Username:
Password:

What's my password?
Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: note [id://1202794]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this?Last hourOther CB clients
Other Users?
Others having a coffee break in the Monastery: (5)
As of 2024-04-23 07:23 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    No recent polls found