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
10 Levenberg-Marquardt fitting is implemented. Other procedures should be
11 added as required. For a fairly concise overview on fitting see
12 Numerical Recipes, chapter 15 "Modeling of data".
13
15 use PDL::Fit::LM;
16 $ym = lmfit $x, $y, $sigma, \&expfunc, $initp, {Maxiter => 300};
17
19 lmfit
20 Levenberg-Marquardt fitting of a user supplied model function
21 ($ym,$finalp,$covar,$iters) =
23 lmfit $x, $y, $sigma, \&expfunc, $initp, {Maxiter => 300, Eps => 1e-3};
24 where $x is the independent variable and $y the value of the dependent
26 variable at each $x, $sigma is the estimate of the uncertainty (i.e.,
27 standard deviation) of $y at each data point, the fourth argument is a
28 subroutine reference (see below), and $initp the initial values of the
29 parameters to be adjusted.
30 Options:
32 Maxiter: maximum number of iterations before giving up
34 Eps: convergence criterion for fit; success when normalized change
35 in chisquare smaller than Eps
36 The user supplied sub routine reference should accept 4 arguments
38 • a vector of independent values $x
40 • a vector of fitting parameters
42 • a vector of dependent variables that will be assigned upon return
44 • a matrix of partial derivatives with respect to the fitting
46 parameters that will be assigned upon return
47 As an example take this definition of a single exponential with 3
49 parameters (width, amplitude, offset):
50 sub expdec {
52 my ($x,$par,$ym,$dyda) = @_;
53 my ($width,$amp,$off) = map {$par->slice("($_)")} (0..2);
54 my $arg = $x/$width;
55 my $ex = exp($arg);
56 $ym .= $amp*$ex+$off;
57 my (@dy) = map {$dyda->slice(",($_)")} (0..2);
58 $dy[0] .= -$amp*$ex*$arg/$width;
59 $dy[1] .= $ex;
60 $dy[2] .= 1;
61 }
62 Note usage of the ".=" operator for assignment
64 In scalar context returns a vector of the fitted dependent variable. In
66 list context returns fitted y-values, vector of fitted parameters, an
67 estimate of the covariance matrix (as an indicator of goodness of fit)
68 and number of iterations performed.
69 An extended example script that uses lmfit is included below. This
71 nice example was provided by John Gehman and should help you to master
72 the initial hurdles. It can also be found in the Example/Fit directory.
73 use PDL;
75 use PDL::Math;
76 use PDL::Fit::LM;
77 use strict;
78 ### fit using pdl's lmfit (Marquardt-Levenberg non-linear least squares fitting)
81 ###
82 ### `lmfit' Syntax:
83 ###
84 ### ($ym,$finalp,$covar,$iters)
85 ### = lmfit $x, $y, $sigma, \&fn, $initp, {Maxiter => 300, Eps => 1e-3};
86 ###
87 ### Explanation of variables
88 ###
89 ### OUTPUT
90 ### $ym = pdl of fitted values
91 ### $finalp = pdl of parameters
92 ### $covar = covariance matrix
93 ### $iters = number of iterations actually used
94 ###
95 ### INPUT
96 ### $x = x data
97 ### $y = y data
98 ### $sigma = piddle of y-uncertainties for each value of $y (can be set to scalar 1 for equal weighting)
99 ### \&fn = reference to function provided by user (more on this below)
100 ### $initp = initial values for floating parameters
101 ### (needs to be explicitly set prior to use of lmfit)
102 ### Maxiter = maximum iterations
103 ### Eps = convergence criterion (maximum normalized change in Chi Sq.)
104 ### Example:
106 # make up experimental data:
107 my $xdata = pdl sequence 5;
108 my $ydata = pdl [1.1,1.9,3.05,4,4.9];
109 # set initial prameters in a pdl (order in accord with fit function below)
111 my $initp = pdl [0,1];
112 # Weight all y data equally (else specify different uncertainties in a pdl)
114 my $sigma = 1;
115 # Use lmfit. Fourth input argument is reference to user-defined
117 # subroutine ( here \&linefit ) detailed below.
118 my ($yf,$pf,$cf,$if) = lmfit $xdata, $ydata, $sigma, \&linefit, $initp;
119 # Note output
121 print "\nXDATA\n$xdata\nY DATA\n$ydata\n\nY DATA FIT\n$yf\n\n";
122 print "Slope and Intercept\n$pf\n\nCOVARIANCE MATRIX\n$cf\n\n";
123 print "NUMBER ITERATIONS\n$if\n\n";
124 # simple example of user defined fit function. Guidelines included on
127 # how to write your own function subroutine.
128 sub linefit {
129 # leave this line as is
131 my ($x,$par,$ym,$dyda) = @_;
132 # $m and $c are fit parameters, internal to this function
134 # call them whatever make sense to you, but replace (0..1)
135 # with (0..x) where x is equal to your number of fit parameters
136 # minus 1
137 my ($m,$c) = map { $par->slice("($_)") } (0..1);
138 # Write function with dependent variable $ym,
140 # independent variable $x, and fit parameters as specified above.
141 # Use the .= (dot equals) assignment operator to express the equality
142 # (not just a plain equals)
143 $ym .= $m * $x + $c;
144 # Edit only the (0..1) part to (0..x) as above
146 my (@dy) = map {$dyda -> slice(",($_)") } (0..1);
147 # Partial derivative of the function with respect to first
149 # fit parameter ($m in this case). Again, note .= assignment
150 # operator (not just "equals")
151 $dy[0] .= $x;
152 # Partial derivative of the function with respect to next
154 # fit parameter ($y in this case)
155 $dy[1] .= 1;
156 # Add $dy[ ] .= () lines as necessary to supply
158 # partial derivatives for all floating parameters.
159 }
160 tlmfit
162 threaded version of Levenberg-Marquardt fitting routine mfit
163 tlmfit $x, $y, float(1)->dummy(0), $na, float(200), float(1e-4),
165 $ym=null, $afit=null, \&expdec;
166 Signature: tlmfit(x(n);y(n);sigma(n);initp(m);iter();eps();[o] ym(n);[o] finalp(m);
168 OtherPar => subref)
169 a threaded version of "lmfit" by using perl threading. Direct threading
171 in "lmfit" seemed difficult since we have an if condition in the
172 iteration. In principle that can be worked around by using "where" but
173 .... Send a threaded "lmfit" version if you work it out!
174 Since we are using perl threading here speed is not really great but it
176 is just convenient to have a threaded version for many applications (no
177 explicit for-loops required, etc). Suffers from some of the current
178 limitations of perl level threading.
179
181 Not known yet.
182
184 This file copyright (C) 1999, Christian Soeller
185 (c.soeller@auckland.ac.nz). All rights reserved. There is no warranty.
186 You are allowed to redistribute this software documentation under
187 certain conditions. For details, see the file COPYING in the PDL
188 distribution. If this file is separated from the PDL distribution, the
189 copyright notice should be included in the file.
190perl v5.32.1 2021-02-15 LM(3)