.
#======================================================================
# . L M F I T . N O R M A L
# doc: Wed Feb 24 09:40:06 1999
# dlm: Fri Jul 28 13:35:24 2006
# (c) 1999 A.M. Thurnherr
# uE-Info: 34 51 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# What you need to provide if you wanna fit a different
# model function to your data:
# - a number of global variables to be set during loading
# - a number of subs to perform admin tasks (usage, init, ...)
# - a sub to evaluate the model function which is to be fitted using
# a number of pararams which are all stored in @A (beginning at
# A[1]!!!). You also need to return the partial derivatives of
# the model function wrt all params.
# - the interface is documented between +++++++ lines
# check if a given distribution is normal
# NB:
# - fitting is based on gauss curve fitting [.lmfit.gauss]
# - heuristics are taken from there and scaled for the normal
# parameter choices
# - simplified, e.g. y-shift is removed (does not make sense for
# distribution)
# - added chi^2 significance testing to &modelCleanup() on -x
# HISTORY:
# Oct 04, 1999: - created from [.lmfit.gauss]
# Oct 05, 1999: - added chi^2 significance test
# - removed -y
# - improved heuristics
# Mar 17, 2001: - param->arg
# Jul 28, 2006: - Version 3.3 [HISTORY]; &isnan()
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
#
# $modelOpts string of allowed options
# $modelOptsUsage usage information string for options
# $modelMinArgs min # of arguments of model
# $modelArgsUsage usage information string for arguments
#
# The following variables may be set later but not after &modelInit()
#
# $modelNFit number of params to fit in model
# @nameA symbolic names of model parameters
#
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
$modelOpts = "x";
$modelOptsUsage = "[-x chi^2 test]";
$modelMinArgs = 0;
$modelArgsUsage = "[area guess [mean guess [sigma guess]]]";
$modelNFit = 3;
$nameA[1] = "area";
$nameA[2] = "mean";
$nameA[3] = "sigma";
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# &modelUsage() mangle parameters; NB: there may be `infinite' # of
# filenames after model arguments; this usually sets
# @A (the model parameters) but these can later be
# calculated heuristically during &modelInit()
#
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
sub modelUsage()
{
$A[1] = nan; $A[2] = nan; $A[3] = nan; # usage
$A[1] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
$A[2] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
$A[3] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
&antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]);
}
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# &modelInit() initializes model after reading of data
#
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
sub modelInit()
{
my($i,$j,$ymin,$ymax,$xatymax);
# --------------------------------------------------
# heuristics for initial model param values
# --------------------------------------------------
$ymin = 1e33, $ymax = -1e33, $xatymax = 0;
for ($i=0; $i<=$#ants_; $i++) {
next if ($antsFlagged[$i]);
$ymin = $ants_[$i][$yfnr]
if ($ants_[$i][$yfnr] < $ymin);
$ymax = $ants_[$i][$yfnr], $xatymax = $ants_[$i][$xfnr]
if ($ants_[$i][$yfnr] > $ymax);
}
$A[1] = $ymax - $ymin if isnan($A[1]); # peak guess
$A[2] = $xatymax if isnan($A[2]); # mean guess
if (isnan($A[3])) { # e-scale guess
for ($i=1;
$i<=$#ants_ && !$antsFlagged[$i]
&& $ants_[$i][$yfnr]-$ymin<0.36*$A[1];
$i++) {}
for ($j=$#ants_;
$j>=1 && !$antsFlagged[$i]
&& $ants_[$j][$yfnr]-$ymin < 0.36*$A[1];
$j--) {}
$A[3] = abs($ants_[$i][$xfnr]-$ants_[$j][$xfnr]) / 2.0;
if ($A[3] == 0.0) {
&antsInfo("$model: sigma heuristic failed (set to 1)!");
$A[3] = 1.0;
}
}
$A[1] *= 1.77 * $A[3] # gauss -> normal
unless (isnan($A[1]) || isnan($A[3]));
$A[3] *= 0.71 unless isnan($A[3]);
}
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# &modelEvaluate(x,A,dyda) evaluate Normal distribution curve at x
# x x value (NOT xfnr)
# A reference to @A
# dyda reference to array for partial derivatives
# (wrt individual params in @A)
# <ret val> y value
#
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
sub modelEvaluate($$$)
{
my($x,$AR,$dydaR) = @_;
my($peak) = $AR->[1] / (2.506628274631 * $AR->[3]);
my($dx ) = $x - $AR->[2];
my($sig2) = $AR->[3] * $AR->[3];
my($expo) = exp(-$dx*$dx/(2*$sig2));
my($norm) = $peak * $expo;
if (defined($dydaR)) {
$dydaR->[1] = $norm / $AR->[1];
$dydaR->[2] = $norm * $dx / $sig2;
$dydaR->[3] = $norm/$AR->[3] * ($dx*$dx/$sig2 - 1);
}
return $norm;
}
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# &modelCleanup() cleans up after fitting but before output
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
sub modelCleanup()
{
return unless ($opt_x);
require "$ANTS/libfuns.pl";
my($chisq) = 0;
my($nval,$prob,$sign);
for ($i=0; $i<=$#ants_; $i++) {
next if ($antsFlagged[$i]);
# next if ($ants_[$i][$yfnr] <= 1); # IGNORE TAIL HEURISTICS
$nval = &modelEvaluate($ants_[$i][$xfnr],\@A);
$chisq += ($ants_[$i][$yfnr] - $nval)**2 / $nval;
}
$prob = &gammq(($ndata-3)/2,$chisq/2);
$sign = int($prob*100);
&antsInfo("$model: normal-distr. hypothesis disproved at $sign%% sign. level");
}