File Coverage

blib/lib/Math/Gauss.pm

Criterion	Covered	Total	%
statement	44	44	100.0
branch	22	24	91.6
condition	3	3	100.0
subroutine	7	7	100.0
pod	0	3	0.0
total	76	81	93.8

package Math::Gauss;

use 5.010000;
use strict;
use warnings;
use Carp;

require Exporter;
our @ISA = qw( Exporter );
our %EXPORT_TAGS = ( 'all' => [qw( pdf cdf inv_cdf )] );
Exporter::export_ok_tags( 'all' );

our $VERSION = '0.01';

my $SQRT2PI = 2.506628274631;

sub pdf {
  my ( $x, $m, $s ) = ( 0, 0, 1 );
  $x = shift if @_;
  $m = shift if @_;
  $s = shift if @_;

  if( $s <= 0 ) {
    croak( "Can't evaluate Math::Gauss:pdf for \$s=$s not strictly positive" );
  }

  my $z = ($x-$m)/$s;

  return exp(-0.5*$z*$z)/($SQRT2PI*$s);
}

sub cdf {
  my ( $x, $m, $s ) = ( 0, 0, 1 );
  $x = shift if @_;
  $m = shift if @_;
  $s = shift if @_;

  # Abramowitz & Stegun, 26.2.17
  # absolute error less than 7.5e-8 for all x

  if( $s <= 0 ) {
    croak( "Can't evaluate Math::Gauss:cdf for \$s=$s not strictly positive" );
  }

  my $z = ($x-$m)/$s;

  my $t = 1.0/(1.0 + 0.2316419*abs($z));
  my $y = $t*(0.319381530
	      + $t*(-0.356563782
		    + $t*(1.781477937
			  + $t*(-1.821255978
				+ $t*1.330274429 ))));
  if( $z > 0 ) {
    return 1.0 - pdf( $z )*$y;
  } else {
    return pdf( $z )*$y;
  }
}

sub inv_cdf {
  my ( $x ) = @_;

  if( $x<=0.0 || $x>=1.0 ) {
    croak( "Can't evaluate Math::Gauss::inv_cdf for \$x=$x outside ]0,1[" );
  }

  # Abramowitz & Stegun, 26.2.23
  # absolute error less than 4.5e-4 for all x

  my $t;
  if( $x < 0.5 ) {
    $t = sqrt( -2.0*log($x) );
  } else {
    $t = sqrt( -2.0*log(1.0 - $x) );
  }

  my $y = (2.515517 + $t*(0.802853 + $t*0.010328));
  $y /= 1.0 + $t*(1.432788 + $t*(0.189269 + $t*0.001308));

  if( $x<0.5 ) {
    return ( $y - $t );
  } else {
    return ( $t - $y );
  }
}

1;
__END__