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'.
This also means we no longer need the normal distribution libraries from lib/

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
 	};
 };
 
diff --git a/LICENSE b/LICENSE
index d17b9c3..ddc491e 100644
--- a/LICENSE
+++ b/LICENSE
@@ -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).
diff --git a/NEWS b/NEWS
index 1d23fd1..23e6969 100644
--- a/NEWS
+++ b/NEWS
@@ -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.