1LM(3) User Contributed Perl Documentation LM(3)
2
6 PDL::Fit::LM -- Levenberg-Marquardt fitting routine for PDL
7
9 This module provides fitting functions for PDL. Currently, only Leven‐
10 berg-Marquardt fitting is implemented. Other procedures should be added
11 as required. For a fairly concise overview on fitting see Numerical
12 Recipes, chapter 15 "Modeling of data".
13
15 use PDL::Fit::LM;
16 $ym = lmfit $x, $y, $sig, \&expfunc, $a, {Maxiter => 300};
17
19 lmfit
20 Levenberg-Marquardt fitting of a user supplied model function
22 ($ym,$a,$covar,$iters) =
24 lmfit $x, $y, $sig, \&expfunc, $a, {Maxiter => 300, Eps => 1e-3};
25 Options:
27 Maxiter: maximum number of iterations before giving up
29 Eps: convergence citerium for fit; success when normalized change
30 in chisquare smaller than Eps
31 The user supplied sub routine reference should accept 4 arguments
33 · a vector of independent values $x
35 · a vector of fitting parameters
37 · a vector of dependent variables that will be assigned upon return
39 · a matrix of partial derivatives with respect to the fitting parame‐
41 ters that will be assigned upon return
42 As an example take this definition of a single exponential with 3
44 parameters (width, amplitude, offset):
45 sub expdec {
47 my ($x,$par,$ym,$dyda) = @_;
48 my ($a,$b,$c) = map {$par->slice("($_)")} (0..2);
49 my $arg = $x/$a;
50 my $ex = exp($arg);
51 $ym .= $b*$ex+$c;
52 my (@dy) = map {$dyda->slice(",($_)")} (0..2);
53 $dy[0] .= -$b*$ex*$arg/$a;
54 $dy[1] .= $ex;
55 $dy[2] .= 1;
56 }
57 Note usage of the ".=" operator for assignment
59 In scalar context returns a vector of the fitted dependent variable. In
61 list context returns fitted y-values, vector of fitted parameters, an
62 estimate of the covariance matrix (as an indicator of goodness of fit)
63 and number of iterations performed.
64 An extended example script that uses lmfit is included below. This
66 nice example was provided by John Gehman and should help you to master
67 the initial hurdles. It can also be found in the Example/Fit directory.
68 use PDL;
70 use PDL::Math;
71 use PDL::Fit::LM;
72 use strict;
73 ### fit using pdl's lmfit (Marquardt-Levenberg non-linear least squares fitting)
75 ###
76 ### `lmfit' Syntax:
77 ###
78 ### ($ym,$a,$covar,$iters)
79 ### = lmfit $x, $y, $sig, \&fn, $initp, {Maxiter => 300, Eps => 1e-3};
80 ###
81 ### Explanation of variables
82 ###
83 ### OUTPUT
84 ### $ym = pdl of fitted values
85 ### $a = pdl of paramters
86 ### $covar = covariance matrix
87 ### $iters = number of iterations actually used
88 ###
89 ### INPUT
90 ### $x = x data
91 ### $y = y data
92 ### $sig = weights for y data (can be set to scalar 1 for equal weighting)
93 ### \&fn = reference to function provided by user (more on this below)
94 ### $initp = initial values for floating parameters
95 ### (needs to be explicitly set prior to use of lmfit)
96 ### Maxiter = maximum iterations
97 ### Eps = convergence criterium (maximum normalized change in Chi Sq.)
98 ### Example:
100 # make up experimental data:
101 my $xdata = pdl sequence 5;
102 my $ydata = pdl [1.1,1.9,3.05,4,4.9];
103 # set initial prameters in a pdl (order in accord with fit function below)
105 my $initp = pdl [0,1];
106 # Weight all y data equally (else specify different weights in a pdl)
108 my $wt = 1;
109 # Use lmfit. Fourth input argument is reference to user-defined
111 # subroutine ( here \&linefit ) detailed below.
112 my ($yf,$pf,$cf,$if) = lmfit $xdata, $ydata, $wt, \&linefit, $initp;
113 # Note output
115 print "\nXDATA\n$xdata\nY DATA\n$ydata\n\nY DATA FIT\n$yf\n\n";
116 print "Slope and Intercept\n$pf\n\nCOVARIANCE MATRIX\n$cf\n\n";
117 print "NUMBER ITERATIONS\n$if\n\n";
118 # simple example of user defined fit function. Guidelines included on
120 # how to write your own function subroutine.
121 sub linefit {
122 # leave this line as is
124 my ($x,$par,$ym,$dyda) = @_;
125 # $m and $b are fit parameters, internal to this function
127 # call them whatever make sense to you, but replace (0..1)
128 # with (0..x) where x is equal to your number of fit parameters
129 # minus 1
130 my ($m,$b) = map { $par->slice("($_)") } (0..1);
131 # Write function with dependent variable $ym,
133 # independent variable $x, and fit parameters as specified above.
134 # Use the .= (dot equals) assignment operator to express the equality
135 # (not just a plain equals)
136 $ym .= $m * $x + $b;
137 # Edit only the (0..1) part to (0..x) as above
139 my (@dy) = map {$dyda -> slice(",($_)") } (0..1);
140 # Partial derivative of the function with respect to first
142 # fit parameter ($m in this case). Again, note .= assignment
143 # operator (not just "equals")
144 $dy[0] .= $x;
145 # Partial derivative of the function with respect to next
147 # fit parameter ($b in this case)
148 $dy[1] .= 1;
149 # Add $dy[ ] .= () lines as necessary to supply
151 # partial derivatives for all floating paramters.
152 }
153 tlmfit
155 threaded version of Levenberg-Marquardt fitting routine mfit
157 tlmfit $x, $y, float(1)->dummy(0), $na, float(200), float(1e-4),
159 $ym=null, $afit=null, \&expdec;
160 Signature:
162 tlmfit(x(n);y(n);sig(n);a(m);iter();eps();[o] ym(n);[o] ao(m);
164 OtherPar => subref)
165 a threaded version of "lmfit" by using perl threading. Direct threading
167 in "lmfit" seemed difficult since we have an if condition in the itera‐
168 tion. In principle that can be worked around by using "where" but ....
169 Send a threaded "lmfit" version if you work it out!
170 Since we are using perl threading here speed is not really great but it
172 is just convenient to have a threaded version for many applications (no
173 explicit for-loops required, etc). Suffers from some of the current
174 limitations of perl level threading.
175
177 Not known yet.
178
180 This file copyright (C) 1999, Christian Soeller (c.soeller@auck‐
181 land.ac.nz). All rights reserved. There is no warranty. You are
182 allowed to redistribute this software documentation under certain con‐
183 ditions. For details, see the file COPYING in the PDL distribution. If
184 this file is separated from the PDL distribution, the copyright notice
185 should be included in the file.
186perl v5.8.8 2006-08-09 LM(3)