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, $sig, \&expfunc, $a, {Maxiter => 300};
17
19 lmfit
20 Levenberg-Marquardt fitting of a user supplied model function
21 ($ym,$a,$covar,$iters) =
23 lmfit $x, $y, $sig, \&expfunc, $a, {Maxiter => 300, Eps => 1e-3};
24 Options:
26 Maxiter: maximum number of iterations before giving up
28 Eps: convergence citerium for fit; success when normalized change
29 in chisquare smaller than Eps
30 The user supplied sub routine reference should accept 4 arguments
32 · a vector of independent values $x
34 · a vector of fitting parameters
36 · a vector of dependent variables that will be assigned upon return
38 · a matrix of partial derivatives with respect to the fitting
40 parameters that will be assigned upon return
41 As an example take this definition of a single exponential with 3
43 parameters (width, amplitude, offset):
44 sub expdec {
46 my ($x,$par,$ym,$dyda) = @_;
47 my ($a,$b,$c) = map {$par->slice("($_)")} (0..2);
48 my $arg = $x/$a;
49 my $ex = exp($arg);
50 $ym .= $b*$ex+$c;
51 my (@dy) = map {$dyda->slice(",($_)")} (0..2);
52 $dy[0] .= -$b*$ex*$arg/$a;
53 $dy[1] .= $ex;
54 $dy[2] .= 1;
55 }
56 Note usage of the ".=" operator for assignment
58 In scalar context returns a vector of the fitted dependent variable. In
60 list context returns fitted y-values, vector of fitted parameters, an
61 estimate of the covariance matrix (as an indicator of goodness of fit)
62 and number of iterations performed.
63 An extended example script that uses lmfit is included below. This
65 nice example was provided by John Gehman and should help you to master
66 the initial hurdles. It can also be found in the Example/Fit directory.
67 use PDL;
69 use PDL::Math;
70 use PDL::Fit::LM;
71 use strict;
72 ### 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
121 # how to write your own function subroutine.
122 sub linefit {
123 # leave this line as is
125 my ($x,$par,$ym,$dyda) = @_;
126 # $m and $b are fit parameters, internal to this function
128 # call them whatever make sense to you, but replace (0..1)
129 # with (0..x) where x is equal to your number of fit parameters
130 # minus 1
131 my ($m,$b) = map { $par->slice("($_)") } (0..1);
132 # Write function with dependent variable $ym,
134 # independent variable $x, and fit parameters as specified above.
135 # Use the .= (dot equals) assignment operator to express the equality
136 # (not just a plain equals)
137 $ym .= $m * $x + $b;
138 # Edit only the (0..1) part to (0..x) as above
140 my (@dy) = map {$dyda -> slice(",($_)") } (0..1);
141 # Partial derivative of the function with respect to first
143 # fit parameter ($m in this case). Again, note .= assignment
144 # operator (not just "equals")
145 $dy[0] .= $x;
146 # Partial derivative of the function with respect to next
148 # fit parameter ($b in this case)
149 $dy[1] .= 1;
150 # Add $dy[ ] .= () lines as necessary to supply
152 # partial derivatives for all floating paramters.
153 }
154 tlmfit
156 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
168 iteration. 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
181 (c.soeller@auckland.ac.nz). All rights reserved. There is no warranty.
182 You are allowed to redistribute this software documentation under
183 certain conditions. For details, see the file COPYING in the PDL
184 distribution. If this file is separated from the PDL distribution, the
185 copyright notice should be included in the file.
186perl v5.12.3 2009-10-17 LM(3)