0
|
1 |
#======================================================================
|
|
2 |
# L F I T . P L
|
|
3 |
# doc: Sat Jul 31 11:24:47 1999
|
20
|
4 |
# dlm: Mon May 18 10:07:47 2015
|
0
|
5 |
# (c) 1999 A.M. Thurnherr
|
20
|
6 |
# uE-Info: 20 61 NIL 0 0 72 2 2 4 NIL ofnI
|
0
|
7 |
#======================================================================
|
|
8 |
|
|
9 |
# LFIT routine from Numerical Recipes adapted to ANTS
|
|
10 |
|
|
11 |
# HISTORY:
|
|
12 |
# Jul 31, 1999: - manually converted from c-source
|
|
13 |
# Aug 01, 1999: - changed funcs() interface
|
|
14 |
# Sep 26, 1999: - made sure right version of covsrt is used
|
|
15 |
# Jun 28, 2001: - re-added commented out code 'cause it's required
|
|
16 |
# on some perl versions
|
|
17 |
# Jan 5, 2012: - BUG: non-numeric x/y were not handled correctly;
|
|
18 |
# this was only easily apparent when the last
|
|
19 |
# record contained non-numeric values
|
20
|
20 |
# May 18, 2015: - added firstRec, lastRec parameters (V6.1)
|
0
|
21 |
|
|
22 |
# Notes:
|
|
23 |
# - x,y,sig are field numbers for data in $ants_
|
|
24 |
# - funcs is passed the current index & xfnr instead of the x-value
|
|
25 |
# - if sig is a negative number, -sig is used as constant input stddev
|
|
26 |
# - @a, @ia, @covar, &funcs passed as refs
|
|
27 |
# - chi square is returned
|
|
28 |
|
|
29 |
require "$ANTS/nrutil.pl";
|
|
30 |
require "$ANTS/covsrt.pl";
|
|
31 |
require "$ANTS/gaussj.pl";
|
|
32 |
|
20
|
33 |
sub lfit(@)
|
0
|
34 |
{
|
20
|
35 |
my($xfnr,$yfnr,$sig,$aR,$iaR,$covarR,$funcsR,$firstRec,$lastRec) = @_;
|
|
36 |
$firstRec = 0 unless defined($firstRec);
|
|
37 |
$lastRec = $#ants_ unless defined($lastRec);
|
0
|
38 |
|
|
39 |
my($i,$j,$k,$l,$m,$mfit); # int
|
|
40 |
my($ym,$wt,$sum,$sig2i,$chisq); # float
|
|
41 |
my(@beta,@afunc); # float[]
|
|
42 |
|
|
43 |
&matrix(\@beta,1,$#{$aR},1,1);
|
|
44 |
&vector(\@afunc,1,$#{$aR});
|
|
45 |
for ($j=1; $j<=$#{$aR}; $j++) {
|
|
46 |
$mfit++ if ($iaR->[$j]);
|
|
47 |
}
|
|
48 |
croak("lfit: no parameters to be fitted") if ($mfit == 0);
|
|
49 |
for ($j=1; $j<=$mfit; $j++) { # REQUIRED FOR SOME PERL VERSIONS!!!
|
|
50 |
for ($k=1;$ k<=$mfit; $k++) {
|
|
51 |
$covarR->[$j][$k] = 0;
|
|
52 |
}
|
|
53 |
$beta[$j][1] = 0;
|
|
54 |
}
|
20
|
55 |
for ($i=$firstRec; $i<=$lastRec; $i++) {
|
0
|
56 |
next if ($antsFlagged[$i]);
|
|
57 |
next unless numberp($ants_[$i][$xfnr]) && numberp($ants_[$i][$yfnr]);
|
|
58 |
&$funcsR($i,$xfnr,\@afunc);
|
|
59 |
$ym = $ants_[$i][$yfnr];
|
|
60 |
if ($mfit < $#{$aR}) {
|
|
61 |
for ($j=1; $j<=$#{$aR}; $j++) {
|
|
62 |
$ym -= $aR->[$j]*$afunc[$j] if (!$iaR->[$j]);
|
|
63 |
}
|
|
64 |
}
|
|
65 |
if ($sig > 0) { # field number
|
|
66 |
$sig2i = 1.0/($ants_[$i][$sig]*$ants_[$i][$sig]);
|
|
67 |
} else { # const value
|
|
68 |
$sig2i = 1.0/($sig*$sig);
|
|
69 |
}
|
|
70 |
for ($j=0,$l=1; $l<=$#{$aR}; $l++) {
|
|
71 |
if ($iaR->[$l]) {
|
|
72 |
$wt = $afunc[$l]*$sig2i;
|
|
73 |
for ($j++,$k=0,$m=1; $m<=$l; $m++) {
|
|
74 |
$covarR->[$j][++$k] += $wt*$afunc[$m] if ($iaR->[$m]);
|
|
75 |
}
|
|
76 |
$beta[$j][1] += $ym*$wt;
|
|
77 |
}
|
|
78 |
}
|
|
79 |
}
|
|
80 |
for ($j=2; $j<=$mfit; $j++) {
|
|
81 |
for ($k=1;$k<$j;$k++) {
|
|
82 |
$covarR->[$k][$j] = $covarR->[$j][$k];
|
|
83 |
# print(STDERR "covarR[$k][$j] = $covarR->[$k][$j]\n");
|
|
84 |
}
|
|
85 |
}
|
|
86 |
&gaussj($covarR,\@beta);
|
|
87 |
for ($j=0,$l=1;$l<=$#{$aR};$l++) {
|
|
88 |
$aR->[$l]=$beta[++$j][1] if ($iaR->[$l]);
|
|
89 |
}
|
20
|
90 |
for ($i=$firstRec; $i<=$lastRec; $i++) {
|
0
|
91 |
next if ($antsFlagged[$i]);
|
|
92 |
next unless numberp($ants_[$i][$xfnr]) && numberp($ants_[$i][$yfnr]);
|
|
93 |
&$funcsR($i,$xfnr,\@afunc);
|
|
94 |
for ($sum=0,$j=1; $j<=$#{$aR}; $j++) {
|
|
95 |
$sum += $aR->[$j]*$afunc[$j];
|
|
96 |
}
|
|
97 |
my($tmpval) = ($sig > 0) ?
|
|
98 |
($ants_[$i][$yfnr] - $sum) / $ants_[$i][$sig] :
|
|
99 |
($ants_[$i][$yfnr] - $sum) / -$sig;
|
|
100 |
$chisq += $tmpval * $tmpval;
|
|
101 |
}
|
|
102 |
&covsrt($covarR,$iaR);
|
|
103 |
return $chisq;
|
|
104 |
}
|
|
105 |
|