diff options
| -rw-r--r-- | Echolot/Stats.pm | 89 | ||||
| -rw-r--r-- | LICENSE | 8 | ||||
| -rw-r--r-- | NEWS | 5 | ||||
| -rw-r--r-- | debian/copyright | 13 | ||||
| -rw-r--r-- | lib/Math/SpecFun/Erf.pm | 565 | ||||
| -rw-r--r-- | lib/Statistics/Distrib/Normal.pm | 438 |
6 files changed, 57 insertions, 1061 deletions
diff --git a/Echolot/Stats.pm b/Echolot/Stats.pm index 905f210..c7255c7 100644 --- a/Echolot/Stats.pm +++ b/Echolot/Stats.pm @@ -1,7 +1,7 @@ package Echolot::Stats; # (c) 2002 Peter Palfrader <peter@palfrader.org> -# $Id: Stats.pm,v 1.56 2003/06/06 09:32:37 weasel Exp $ +# $Id: Stats.pm,v 1.57 2003/06/06 11:27:59 weasel Exp $ # =pod @@ -21,12 +21,6 @@ use strict; use English; use Echolot::Log; -use Statistics::Distrib::Normal qw{}; - -my $NORMAL = new Statistics::Distrib::Normal; -$NORMAL->mu(0); -$NORMAL->sigma(1); - my $STATS_DAYS; my $SECONDS_PER_DAY; @@ -190,11 +184,36 @@ sub build_list2_capsstr($) { return $str; } +sub median($) { + my ($arr) = @_; + my $cnt = scalar @$arr; + if ($cnt == 0) { + return undef; + } elsif ($cnt % 2 == 0) { + return ($arr->[ int(($cnt - 1 ) / 2) ] + $arr->[ int($cnt / 2) ] ) / 2; + } else { + return $arr->[ int(($cnt - 1 ) / 2) ]; + }; +}; + +sub percentile($$) { + my ($lat, $lats) = @_; + + my $num = scalar @$lats; + my $i; + for (my $i=0; $i < $num; $i++) { + last if $lat < $lats->[$i]; + } + return ($num - $i) / $num; +} sub calculate($$) { my ($addr, $types) = @_; my $now = time(); + my $SKEW_ABS = 15*60; + my $SKEW_PERCENT = 0.80; + my @out; my @done; @@ -207,51 +226,47 @@ sub calculate($$) { }; }; - my $latency = 0; - my $received = 0; - my $sent = 0; - my @latency; - my @received; - my @sent; + my @latency_total = map { $_->[1] } @done; + my @latency_day; + my $sent_total; + my $received_total; + my @sent_day; + my @received_day; for my $done (@done) { - $latency += $done->[1]; $latency [int(($now - $done->[0]) / $SECONDS_PER_DAY)] += $done->[1]; - $sent ++; $sent [int(($now - $done->[0]) / $SECONDS_PER_DAY)] ++; - $received ++; $received[int(($now - $done->[0]) / $SECONDS_PER_DAY)] ++; - }; - $latency /= (scalar @done) if (scalar @done); - $latency = undef unless (scalar @done); - for ( 0 .. $STATS_DAYS - 1 ) { - $latency[$_] /= $received[$_] if ($received[$_]); + push @{ $latency_day [int(($now - $done->[0]) / $SECONDS_PER_DAY)] }, $done->[1]; + $sent_total ++; $sent_day [int(($now - $done->[0]) / $SECONDS_PER_DAY)] ++; + $received_total ++; $received_day[int(($now - $done->[0]) / $SECONDS_PER_DAY)] ++; }; - my $variance = 0; - $variance += ($latency - $_->[1]) ** 2 for (@done); - $variance /= (scalar @done) if (scalar @done); - - my $deviation = sqrt($variance); + @latency_total = sort { $a <=> $b } @latency_total; + my $latency_median = median (\@latency_total); + my @latency_median_day; + for ( 0 .. $STATS_DAYS - 1 ) { + @{$latency_day[$_]} = sort { $a <=> $b } $latency_day[$_]; + $latency_median_day[$_] = median ( $latency_day[$_] ); + } if (scalar @out) { - my @p = - ($deviation != 0) ? - $NORMAL->utp( map { ($now - $_ - $latency) / $deviation } @out ) : + my @p = ( scalar @latency_total ) ? + map { percentile( ($now - $_ - $SKEW_ABS)/$SKEW_PERCENT , \@latency_total ) } @out : map { 0 } @out; for (my $i=0; $i < scalar @out; $i++) { - $sent ++; $sent [int(($now - $out[$i]) / $SECONDS_PER_DAY)] ++; - $received += $p[$i]; $received[int(($now - $out[$i]) / $SECONDS_PER_DAY)] += $p[$i]; + $sent_total ++; $sent_day [int(($now - $out[$i]) / $SECONDS_PER_DAY)] ++; + $received_total += $p[$i]; $received_day[int(($now - $out[$i]) / $SECONDS_PER_DAY)] += $p[$i]; }; }; - $received /= $sent if ($sent); + $received_total /= $sent_total if ($sent_total); for ( 0 .. $STATS_DAYS - 1 ) { - $received[$_] /= $sent[$_] if ($sent[$_]); + $received_day[$_] /= $sent_day[$_] if ($sent_day[$_]); }; return { - avr_latency => $latency, - avr_reliability => $received, - latency_day => \@latency, - reliability_day => \@received + avr_latency => $latency_median, + avr_reliability => $received_total, + latency_day => \@latency_day, + reliability_day => \@received_day }; }; @@ -13,11 +13,3 @@ GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - - -The lib directory includes software by Peter J. Acklam <pjacklam@online.no>: - - Math::SpecFun::Erf - - Statistics::Distrib::Normal - -Those modules may be distributed under the same terms as Perl itself (the GNU -General Public License and the Artistic License). @@ -13,6 +13,11 @@ Changes in shut down in stop and restart. * Regularily clean up the temp directory. * Send out pings every two hours by default. + * Modified statistics generation: + We no longer assume a normal distribution of latencies but instead + use percentiles when calculating a life probability of an outstanding + ping. Also we do not show the mean of latency but the median as + this seems to be 'more correct'. Changes in version * Minor documentation fixes suggested by Ryan Lackey. diff --git a/debian/copyright b/debian/copyright index ed455ab..971dea6 100644 --- a/debian/copyright +++ b/debian/copyright @@ -23,16 +23,3 @@ Place, Suite 330, Boston, MA 02111-1307 USA On Debian systems you can find a copy of the GNU General Public License in /usr/share/common-licenses/GPL. - - - -This package (and the echolot upstream tarball) also includes software by -Peter J. Acklam <pjacklam@online.no>: - - - Math::SpecFun::Erf - - Statistics::Distrib::Normal - -Those modules may be distributed under the same terms as Perl itself (the GNU -General Public License and the Artistic License). - -You can find a copy of the Artistic License in /usr/share/common-licenses too. diff --git a/lib/Math/SpecFun/Erf.pm b/lib/Math/SpecFun/Erf.pm deleted file mode 100644 index c913792..0000000 --- a/lib/Math/SpecFun/Erf.pm +++ /dev/null @@ -1,565 +0,0 @@ -# -# Author: Peter J. Acklam -# Time-stamp: 2000-11-29 23:04:53 -# E-mail: pjacklam@online.no -# URL: http://home.online.no/~pjacklam - -=head1 NAME - -Math::SpecFun::Erf - error and scaled and unscaled complementary error -functions and their inverses - -=head1 SYNOPSIS - - use Math::SpecFun::Erf qw(erf erfc erfcx erfinv erfcinv erfcxinv); - -Imports all the routines explicitly. Use a subset of the list for the -routines you want. - - use Math::SpecFun::Erf qw(:all); - -Imports all the routines, as well. - -=head1 DESCRIPTION - -This module implements the error function, C<erf>, and its inverse -C<erfinv>, the complementary error function, C<erfc>, and its inverse -C<erfcinv>, and the scaled complementary error function, C<erfcx>, and its -inverse C<erfcxinv>. - -For references and details about the algorithms, see the comments inside -this module. - -=head1 FUNCTIONS - -=over 8 - -=item erf EXPR - -=item erf - -Returns the error function evaluated at EXPR. If EXPR is omitted, C<$_> is -used. The error function is - - erf(x) = 2/sqrt(PI) * integral from 0 to x of exp(-t*t) dt - -=item erfinv EXPR - -=item erfinv - -Returns the inverse of the error function evaluated at EXPR. If EXPR is -omitted, C<$_> is used. - -=item erfc EXPR - -=item erfc - -Returns the complementary error function evaluated at EXPR. If EXPR is -omitted, C<$_> is used. The complementary error function is - - erfc(x) = 2/sqrt(PI) * integral from x to infinity of exp(-t*t) dt - = 1 - erf(x) - -Here is a function returning the lower tail probability of the standard -normal distribution function - - use Math::SpecFun::Erf qw(erfc); - - sub ltpnorm ($) { - erfc( - $_[0] / sqrt(2) )/2; - } - -=item erfcinv EXPR - -=item erfcinv - -Returns the inverse complementary error function evaluated at EXPR. If EXPR -is omitted, C<$_> is used. - -Here is a function returning the lower tail quantile of the standard normal -distribution function - - use Math::SpecFun::Erf qw(erfcinv); - - sub ltqnorm ($) { - -sqrt(2) * erfcinv( 2 * $_[0] ); - } - -=item erfcx EXPR - -=item erfcx - -Returns the scaled complementary error function evaluated at EXPR. If EXPR -is omitted, C<$_> is used. The scaled complementary error function is - - erfcx(x) = exp(x*x) * erfc(x) - -=item erfcxinv EXPR - -=item erfcxinv - -Returns the inverse scaled complementary error function evaluated at EXPR. -If EXPR is omitted, C<$_> is used. - -=back - -=head1 HISTORY - -=over 4 - -=item Version 0.03 - -Added the inverse functions. - -=item Version 0.02 - -Minor code tweaking. - -=item Version 0.01 - -First release. - -=back - -=head1 AUTHOR - -Perl translation by Peter J. Acklam E<lt>pjacklam@online.noE<gt> - -FORTRAN code by W. J. Cody, Argonne National Laboratory, March 19, 1990. -FORTRAN code can be found at http://www.netlib.org/specfun/erf - -=head1 COPYRIGHT - -Copyright (c) 1999-2000 Peter J. Acklam. All rights reserved. -This program is free software; you can redistribute it and/or -modify it under the same terms as Perl itself. - -=cut - -package Math::SpecFun::Erf; -require 5.000; -require Exporter; - -use strict; -use vars qw($VERSION @ISA @EXPORT_OK %EXPORT_TAGS); - -$VERSION = '0.02'; -@ISA = qw(Exporter); -@EXPORT_OK = qw(erf erfc erfcx erfinv erfcinv erfcxinv); -%EXPORT_TAGS = ( all => [ @EXPORT_OK ] ); - -######################################################################## -## Internal functions. -######################################################################## - -sub calerf { - my ($arg, $result, $jint) = @_; - local $[ = 1; -#------------------------------------------------------------------ -# -# This packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x) -# for a real argument x. It contains three FUNCTION type -# subprograms: ERF, ERFC, and ERFCX (or DERF, DERFC, and DERFCX), -# and one SUBROUTINE type subprogram, CALERF. The calling -# statements for the primary entries are: -# -# Y=ERF(X) (or Y=DERF(X)), -# -# Y=ERFC(X) (or Y=DERFC(X)), -# and -# Y=ERFCX(X) (or Y=DERFCX(X)). -# -# The routine CALERF is intended for internal packet use only, -# all computations within the packet being concentrated in this -# routine. The function subprograms invoke CALERF with the -# statement -# -# CALL CALERF(ARG,RESULT,JINT) -# -# where the parameter usage is as follows -# -# Function Parameters for CALERF -# call ARG Result JINT -# -# ERF(ARG) ANY REAL ARGUMENT ERF(ARG) 0 -# ERFC(ARG) ABS(ARG) < XBIG ERFC(ARG) 1 -# ERFCX(ARG) XNEG < ARG < XMAX ERFCX(ARG) 2 -# -# The main computation evaluates near-minimax approximations -# from "Rational Chebyshev approximations for the error function" -# by W. J. Cody, Math. Comp., 1969, PP. 631-638. This -# transportable program uses rational functions that theoretically -# approximate erf(x) and erfc(x) to at least 18 significant -# decimal digits. The accuracy achieved depends on the arithmetic -# system, the compiler, the intrinsic functions, and proper -# selection of the machine-dependent constants. -# -#******************************************************************* -#******************************************************************* -# -# Explanation of machine-dependent constants -# -# XMIN = the smallest positive floating-point number. -# XINF = the largest positive finite floating-point number. -# XNEG = the largest negative argument acceptable to ERFCX; -# the negative of the solution to the equation -# 2*exp(x*x) = XINF. -# XSMALL = argument below which erf(x) may be represented by -# 2*x/sqrt(pi) and above which x*x will not underflow. -# A conservative value is the largest machine number X -# such that 1.0 + X = 1.0 to machine precision. -# XBIG = largest argument acceptable to ERFC; solution to -# the equation: W(x) * (1-0.5/x**2) = XMIN, where -# W(x) = exp(-x*x)/[x*sqrt(pi)]. -# XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to -# machine precision. A conservative value is -# 1/[2*sqrt(XSMALL)] -# XMAX = largest acceptable argument to ERFCX; the minimum -# of XINF and 1/[sqrt(pi)*XMIN]. -# -# Approximate values for some important machines are: -# -# XMIN XINF XNEG XSMALL -# -# CDC 7600 (S.P.) 3.13E-294 1.26E+322 -27.220 7.11E-15 -# CRAY-1 (S.P.) 4.58E-2467 5.45E+2465 -75.345 7.11E-15 -# IEEE (IBM/XT, -# SUN, etc.) (S.P.) 1.18E-38 3.40E+38 -9.382 5.96E-8 -# IEEE (IBM/XT, -# SUN, etc.) (D.P.) 2.23D-308 1.79D+308 -26.628 1.11D-16 -# IBM 195 (D.P.) 5.40D-79 7.23E+75 -13.190 1.39D-17 -# UNIVAC 1108 (D.P.) 2.78D-309 8.98D+307 -26.615 1.73D-18 -# VAX D-Format (D.P.) 2.94D-39 1.70D+38 -9.345 1.39D-17 -# VAX G-Format (D.P.) 5.56D-309 8.98D+307 -26.615 1.11D-16 -# -# -# XBIG XHUGE XMAX -# -# CDC 7600 (S.P.) 25.922 8.39E+6 1.80X+293 -# CRAY-1 (S.P.) 75.326 8.39E+6 5.45E+2465 -# IEEE (IBM/XT, -# SUN, etc.) (S.P.) 9.194 2.90E+3 4.79E+37 -# IEEE (IBM/XT, -# SUN, etc.) (D.P.) 26.543 6.71D+7 2.53D+307 -# IBM 195 (D.P.) 13.306 1.90D+8 7.23E+75 -# UNIVAC 1108 (D.P.) 26.582 5.37D+8 8.98D+307 -# VAX D-Format (D.P.) 9.269 1.90D+8 1.70D+38 -# VAX G-Format (D.P.) 26.569 6.71D+7 8.98D+307 -# -#******************************************************************* -#******************************************************************* -# -# Error returns -# -# The program returns ERFC = 0 for ARG >= XBIG; -# -# ERFCX = XINF for ARG < XNEG; -# and -# ERFCX = 0 for ARG >= XMAX. -# -# -# Intrinsic functions required are: -# -# ABS, AINT, EXP -# -# -# Author: W. J. Cody -# Mathematics and Computer Science Division -# Argonne National Laboratory -# Argonne, IL 60439 -# -# Latest modification: March 19, 1990 -# -# Translation to Perl by Peter J. Acklam, December 3, 1999 -# -#------------------------------------------------------------------ - my ($i); - my ($x, $del, $xden, $xnum, $y, $ysq); -#------------------------------------------------------------------ -# Mathematical constants -#------------------------------------------------------------------ - my ($four, $one, $half, $two, $zero) = (4, 1, 0.5, 2, 0); - my $sqrpi = 5.6418958354775628695e-1; - my $thresh = 0.46875; - my $sixten = 16; -#------------------------------------------------------------------ -# Machine-dependent constants -#------------------------------------------------------------------ - my ($xinf, $xneg, $xsmall) = (1.79e308, -26.628, 1.11e-16); - my ($xbig, $xhuge, $xmax) = (26.543, 6.71e7, 2.53e307); -#------------------------------------------------------------------ -# Coefficients for approximation to erf in first interval -#------------------------------------------------------------------ - my @a = (3.16112374387056560e00, 1.13864154151050156e02, - 3.77485237685302021e02, 3.20937758913846947e03, - 1.85777706184603153e-1); - my @b = (2.36012909523441209e01, 2.44024637934444173e02, - 1.28261652607737228e03, 2.84423683343917062e03); -#------------------------------------------------------------------ -# Coefficients for approximation to erfc in second interval -#------------------------------------------------------------------ - my @c = (5.64188496988670089e-1, 8.88314979438837594e0, - 6.61191906371416295e01, 2.98635138197400131e02, - 8.81952221241769090e02, 1.71204761263407058e03, - 2.05107837782607147e03, 1.23033935479799725e03, - 2.15311535474403846e-8); - my @d = (1.57449261107098347e01, 1.17693950891312499e02, - 5.37181101862009858e02, 1.62138957456669019e03, - 3.29079923573345963e03, 4.36261909014324716e03, - 3.43936767414372164e03, 1.23033935480374942e03); -#------------------------------------------------------------------ -# Coefficients for approximation to erfc in third interval -#------------------------------------------------------------------ - my @p = (3.05326634961232344e-1, 3.60344899949804439e-1, - 1.25781726111229246e-1, 1.60837851487422766e-2, - 6.58749161529837803e-4, 1.63153871373020978e-2); - my @q = (2.56852019228982242e00, 1.87295284992346047e00, - 5.27905102951428412e-1, 6.05183413124413191e-2, - 2.33520497626869185e-3); -#------------------------------------------------------------------ - $x = $arg; - $y = abs($x); - if ($y <= $thresh) { -#------------------------------------------------------------------ -# Evaluate erf for |X| <= 0.46875 -#------------------------------------------------------------------ - $ysq = $zero; - if ($y > $xsmall) { $ysq = $y * $y } - $xnum = $a[5]*$ysq; - $xden = $ysq; - for (my $i = 1 ; $i <= 3 ; ++$i) { - $xnum = ($xnum + $a[$i]) * $ysq; - $xden = ($xden + $b[$i]) * $ysq; - } - $$result = $x * ($xnum + $a[4]) / ($xden + $b[4]); - if ($jint != 0) { $$result = $one - $$result } - if ($jint == 2) { $$result = exp($ysq) * $$result } - goto x800; -#------------------------------------------------------------------ -# Evaluate erfc for 0.46875 <= |X| <= 4.0 -#------------------------------------------------------------------ - } elsif ($y <= $four) { - $xnum = $c[9]*$y; - $xden = $y; - for (my $i = 1 ; $i <= 7 ; ++$i) { - $xnum = ($xnum + $c[$i]) * $y; - $xden = ($xden + $d[$i]) * $y; - } - $$result = ($xnum + $c[8]) / ($xden + $d[8]); - if ($jint != 2) { - $ysq = int($y*$sixten)/$sixten; - $del = ($y-$ysq)*($y+$ysq); - $$result = exp(-$ysq*$ysq) * exp(-$del) * $$result; - } -#------------------------------------------------------------------ -# Evaluate erfc for |X| > 4.0 -#------------------------------------------------------------------ - } else { - $$result = $zero; - if ($y >= $xbig) { - if (($jint != 2) || ($y >= $xmax)) { goto x300 } - if ($y >= $xhuge) { - $$result = $sqrpi / $y; - goto x300; - } - } - $ysq = $one / ($y * $y); - $xnum = $p[6]*$ysq; - $xden = $ysq; - for (my $i = 1 ; $i <= 4 ; ++$i) { - $xnum = ($xnum + $p[$i]) * $ysq; - $xden = ($xden + $q[$i]) * $ysq; - } - $$result = $ysq *($xnum + $p[5]) / ($xden + $q[5]); - $$result = ($sqrpi - $$result) / $y; - if ($jint != 2) { - $ysq = int($y*$sixten)/$sixten; - $del = ($y-$ysq)*($y+$ysq); - $$result = exp(-$ysq*$ysq) * exp(-$del) * $$result; - } - } -#------------------------------------------------------------------ -# Fix up for negative argument, erf, etc. -#------------------------------------------------------------------ - x300: - if ($jint == 0) { - $$result = ($half - $$result) + $half; - if ($x < $zero) { $$result = -$$result } - } elsif ($jint == 1) { - if ($x < $zero) { $$result = $two - $$result } - } else { - if ($x < $zero) { - if ($x < $xneg) { - $$result = $xinf; - } else { - $ysq = int($x*$sixten)/$sixten; - $del = ($x-$ysq)*($x+$ysq); - $y = exp($ysq*$ysq) * exp($del); - $$result = ($y+$y) - $$result; - } - } - } - x800: - return 1; -#---------- Last card of CALERF ---------- -} - -sub erf { - my $x = @_ ? $_[0] : $_; -#-------------------------------------------------------------------- -# -# This subprogram computes approximate values for erf(x). -# (see comments heading CALERF). -# -# Author/date: W. J. Cody, January 8, 1985 -# -# Translation to Perl by Peter J. Acklam, December 3, 1999 -# -#-------------------------------------------------------------------- - my $result; - my $jint = 0; - calerf($x, \$result, $jint); - my $erf = $result; - return $erf; -#---------- Last card of ERF ---------- -} - -######################################################################## -## User functions. -######################################################################## - -sub erfc { - my $x = @_ ? $_[0] : $_; -#-------------------------------------------------------------------- -# -# This subprogram computes approximate values for erfc(x). -# (see comments heading CALERF). -# -# Author/date: W. J. Cody, January 8, 1985 -# -# Translation to Perl by Peter J. Acklam, December 3, 1999 -# -#-------------------------------------------------------------------- - my ($result); - my $jint = 1; - calerf($x, \$result, $jint); - my $erfc = $result; - return $erfc; -#---------- Last card of ERFC ---------- -} - -sub erfcx { - my $x = @_ ? $_[0] : $_; -#------------------------------------------------------------------ -# -# This subprogram computes approximate values for exp(x*x) * erfc(x). -# (see comments heading CALERF). -# -# Author/date: W. J. Cody, March 30, 1987 -# -# Translation to Perl by Peter J. Acklam, December 3, 1999 -# -#------------------------------------------------------------------ - my ($result); - my $jint = 2; - calerf($x, \$result, $jint); - my $erfcx = $result; - return $erfcx; -#---------- Last card of ERFCX ---------- -} - -sub erfinv { - my $y = @_ ? $_[0] : $_; - - return 0 if $y == 0; - return erfcinv(1-$y) if $y > 0.5; - return -erfcinv(1+$y) if $y < -0.5; - - # - # Halley's rational 3rd order method: - # u <- f(x)/f'(x) - # v <- f''(x)/f'(x) - # x <- x - u/(1-u*v/2) - # - # Here: - # f(x) = erf(x) - y - # f'(x) = 2/sqrt(pi)*exp(-x*x) - # f''(x) = -4/sqrt(pi)*x*exp(-x*x) - # - my $x = 0; - my $dx; - my $c = .88622692545275801364908374167055; # sqrt(pi)/2 - my $eps = 5e-15; - do { - my $f = erf($x) - $y; - my $u = $c*$f*exp($x*$x); - $dx = -$u/(1+$u*$x); - $x += $dx; - } until abs($dx/$x) <= $eps; - return $x; -} - -sub erfcinv { - my $y = @_ ? $_[0] : $_; - - return 0 if $y == 1; - - my $flag = $y > 1; - $y = 2 - $y if $flag; - - # - # Halley's rational 3rd order method: - # u <- f(x)/f'(x) - # v <- f''(x)/f'(x) - # x <- x - u/(1-u*v/2) - # - # Here: - # f(x) = erfc(x) - y - # f'(x) = -2/sqrt(pi)*exp(-x*x) - # f''(x) = 4/sqrt(pi)*x*exp(-x*x) - # - my $x = 0; - my $dx; - my $c = -.88622692545275801364908374167055; # sqrt(pi)/2 - my $eps = 5e-15; - do { - my $u = $c*(erfcx($x) - $y*exp($x*$x)); - $dx = -$u/(1+$u*$x); - $x += $dx; - } until abs($dx/$x) <= $eps; - - return $flag ? -$x : $x; -} - -sub erfcxinv { - my $y = @_ ? $_[0] : $_; - - return 0 if $y == 1; - - # - # Halley's 3rd order method: - # u <- f(x)/f'(x) - # v <- f''(x)/f'(x) - # x <- x - u/(1-u*v/2) - # - # Here: - # f(x) = erfcx(x) - y - # f'(x) = 2*(x*erfcx(x)-1/sqrt(pi)); - # f''(x) = (2+4*x*x)*erfcx(x) - 4*x/sqrt(pi); - # - my $x = 0; - my $dx; - my $c = .56418958354775628694807945156079; # 1/sqrt(pi) - my $d = 2.2567583341910251477923178062432; # 4/sqrt(pi) - my $eps = 5e-15; - do { - my $f = erfcx($x) - $y; - my $df = 2*($x*erfcx($x)-$c); - my $ddf = (2+4*$x*$x)*erfcx($x) - $x*$d; - my $u = $f/$df; - my $v = $ddf/$df; - $dx = -$u/(1-$u*$v/2); - $x += $dx; - } until abs($dx/$x) <= $eps; - return $x; -} diff --git a/lib/Statistics/Distrib/Normal.pm b/lib/Statistics/Distrib/Normal.pm deleted file mode 100644 index 58cf737..0000000 --- a/lib/Statistics/Distrib/Normal.pm +++ /dev/null @@ -1,438 +0,0 @@ -# -# Author: Peter J. Acklam -# Time-stamp: 2000-06-02 23:42:57 -# E-mail: pjacklam@online.no -# URL: http://home.online.no/~pjacklam - -=head1 NAME - -Statistics::Distrib::Normal - the normal distribution - -=head1 SYNOPSIS - - use Statistics::Distrib::Normal; - - $dist = new Statistics::Distrib::Normal; - - $dist->mu(3); # set the location parameter - $dist->sigma(5); # set the scale parameter - @x = $dist->rand(10); # generate random numbers - - # or - - @x = Statistics::Distrib::Normal->new(Mu => 3, Sigma => 5)->rand(10); - -=head1 DESCRIPTION - -This module contains miscellaneous functions related to the normal -distribution. - -=cut - -package Statistics::Distrib::Normal; -require 5.000; - -use strict; -use vars qw($VERSION); -use Carp; - -$VERSION = '0.01'; - -use constant PI => 4 * atan2 1, 1; -use constant TWOPI => 2 * PI; - -# the smallest positive floating-point number such that 1+EPS > 1 -use constant EPS => 2.220446049250313080847263336181640625e-016; - -## -## Constructor -## - -=head1 CONSTRUCTOR - -=over 4 - -=item new ( [ OPTIONS ] ) - -C<OPTIONS> is a list of options given in the form of key-value -pairs, just like a hash table. Valid options are - -=over 8 - -=item B<Mu> - -Sets the mu parameter (the mean) of the distribution to the specified -value. - -=item B<Sigma> - -Sets the sigma parameter (the standard deviation) of the distribution -to the specified value. - -=back - -=back - -=cut - -sub new { - my $self = shift; - my $class = ref($self) || $self; - my %arg = @_; - - my %hash = ( mu => 0, - sigma => 1, - ); - - my $me = bless \%hash, $class; - - foreach my $key ( keys %arg ) { - $me->mu($arg{Mu}), next if $key eq 'Mu'; - $me->sigma($arg{Sigma}), next if $key eq 'Sigma'; - carp "Unknown option $arg{$key} ignored"; - } - - return $me; -} - -## -## Methods -## - -=pod - -=head1 METHODS - -=over 4 - -=item mean ( [ MEAN ] ) - -Set the mu parameter (the mean) of the distribution to C<MU>. If -C<MU> is omitted, the current value of mu is returned. - -=cut - -sub mu { - my $me = shift; - croak 'Too many arguments' if @_ > 1; - if ( @_ ) { - $me->{mu} = shift; - return 1; - } - return $me->{mu}; -} - -=pod - -=item sigma ( [ SDEV ] ) - -Set the sigma parameter (the standard deviation) of the distribution -to C<SIGMA>. If C<SIGMA> is omitted, the current value of sigma is -returned. - -=cut - -sub sigma { - my $me = shift; - croak 'Too many arguments' if @_ > 1; - if ( @_ ) { - my $sigma = shift; - croak 'Standard deviation most be positive' unless $sigma > 0; - $me->{sigma} = $sigma; - return 1; - } - return $me->{sigma}; -} - -=pod - -=item pdf ( X1 [, X2 [, X3 ... ] ] ) - -Evaluate the probability density function at C<X1>, C<X2>, C<X3>, ... - -=cut - -sub pdf { - my $me = shift; - croak 'Not enough arguments' unless @_; - my $mu = $me->{mu}; - my $sigma = $me->{sigma}; - my $const = log(TWOPI * $sigma * $sigma); - my @f; - foreach my $x ( @_ ) { - my $z = ( $x - $mu ) / $sigma; - push @f, exp( -0.5 * ( $const + $z*$z ) ); - } - return @f; -} - -=pod - -=item ltp ( X1 [, X2 [, ... ] ] ) - -Evaluate the lower tail probability function at C<X1>, C<X2>, C<X3>, -... - -=cut - -sub ltp { - my $me = shift; - croak 'Not enough arguments' unless @_; - my $mu = $me->{mu}; - my $sigma = $me->{sigma}; - - require Math::SpecFun::Erf; - import Math::SpecFun::Erf qw(erfc); - - my @p; - foreach my $x ( @_ ) { - my $z = ( $x - $mu ) / $sigma; - push @p, erfc( - $_[0] / sqrt(2) )/2; - } - return @p; -} - -=pod - -=item utp ( X1 [, X2 [, ... ] ] ) - -Evaluate the upper tail probability function at C<X1>, C<X2>, C<X3>, -... - -=cut - -sub utp { - my $me = shift; - croak 'Not enough arguments' unless @_; - my $mu = $me->{mu}; - my $sigma = $me->{sigma}; - - require Math::SpecFun::Erf; - import Math::SpecFun::Erf qw(erfc); - - my @p; - foreach my $x ( @_ ) { - my $z = ( $x - $mu ) / $sigma; - push @p, erfc( $_[0] / sqrt(2) )/2; - } - return @p; -} - -=pod - -=item ltq ( P1 [, P2 [, ... ] ] ) - -Returns the lower tail quantile for the probabilities C<P1>, C<P2>, -C<P3>, ... - -=cut - -sub ltq { - croak 'Method not implemented yet'; -} - -=pod - -=item utq ( P1 [, P2 [, P3 ... ] ] ) - -Returns the upper tail quantile for the probabilities C<P1>, C<P2>, -C<P3>, ... - -=cut - -sub utq { - croak 'Method not implemented yet'; -} - -=pod - -=item intprob( XLOW, XHIGH ) - -Interval probability. Returns the probability of an outcome in the -interval from XLOW to XHIGH. - -=cut - -sub intprob { - my $me = shift; - croak 'Bad number of arguments' unless @_ == 2; - my ($xlow, $xhigh) = @_; - return 0 unless $xlow < $xhigh; - my $mu = $me->{mu}; - my $sigma = $me->{sigma}; - - if ( $mu < $xlow ) { - return $me->utp($xlow) - $me->utp($xhigh); - } else { - return $me->ltp($xhigh) - $me->ltp($xlow); - } -} - -=pod - -=item rand( [ NUM ] ) - -Generate random variables. If C<NUM> is omitted, a single variable is -returned. - -=cut - -sub rand { - my $me = shift; - my $num; - if ( @_ ) { - $num = shift; - croak 'Too many arguments' if @_; - croak 'Argument must be positive integer' - unless ($num == int $num) && ($num > 0); - } else { - $num = 1; - } - - # Generate the random variables by the Box-Muller method. - my @z; - my $mu = $me->{mu}; - my $sigma = $me->{sigma}; - my $const = -2 * $sigma * $sigma; - my $i; - for ( $i = 0 ; $i < $num ; $i += 2 ) { - my $r = sqrt $const * log rand; - my $t = TWOPI * rand; - push @z, $mu + $r * sin $t, $mu + $r * cos $t; - } - pop @z if $i > $num; - return @z; -} - -=pod - -=item expectation () - -Return the expectation of the distribution. - -=cut - -sub expectation { - my $me = shift; - croak 'Too many arguments' if @_; - return $me->{mu}; -} - -=pod - -=item variance () - -Return the variance of the distribution. - -=cut - -sub variance { - my $me = shift; - croak 'Too many arguments' if @_; - return $me->{sigma}**2; -} - -=pod - -=item skewness () - -Return the skewness of the distribution. - -=cut - -sub skewness { - my $me = shift; - croak 'Too many arguments' if @_; - return 0; -} - -=pod - -=item kurtosis () - -Return the kurtosis (normalized) of the distribution. - -=cut - -sub kurtosis { - my $me = shift; - croak 'Too many arguments' if @_; - return 0; -} - -=item dmo - -Direct moments for the distribution. - -Not implemented yet. - -=cut - -sub dmo { - croak 'Method not implemented yet'; -} - -=pod - -=item cmo - -Central moments for the distribution. - -=cut - -sub cmo { - croak 'Method not implemented yet'; -} - -=pod - -=item mode () - -Returns the mode of the distribution. - -=cut - -sub mode { - my $me = shift; - croak 'Too many input arguments' if @_; - return $me->{mu}; -} - -=back - -=head1 BUGS - -None known. - -=head1 LIMITATIONS - -Degenerate cases are not allowed for most methods; e.g., a -distribution with zero variance. - -=head1 HISTORY - -=over 4 - -=item Version 0.02 - -Code formatting changes. - -=item Version 0.01 - -First release. - -=back - -=head1 AUTHOR - -Peter J. Acklam E<lt>pjacklam@online.noE<gt>. - -=head1 COPYRIGHT/LICENSE - -Copyright (c) 1999-2000 Peter J. Acklam. All rights reserved. -This program is free software; you can redistribute it and/or -modify it under the same terms as Perl itself. - -=cut - -1; # Modules must return true. |