09/11/27 11:57:05
my $probability = normdist( x => 160, mean => 170, sigma => 10, precision => 10 );
print $probability;
sub normdist {
my %param = @_;
use Math::Trig qw( pi );
my $h = 1/1_000;
my $precision = 10**-$param{precision};
my $z = ( $param{x} - $param{mean} ) / $param{sigma};
my $const = 1 / sqrt( 2*pi() );
my $fx = sub { $const * exp( -0.5 * $_[0]**2 ) };
my $s = 0;
my $n = 0;
while (1) {
my $y1 = $fx->( $z + $n*$h );
my $y2 = $fx->( $z + ($n-1)*$h );
my $y3 = $fx->( $z + ($n-2)*$h );
my $ds = $h/3 * ($y1 + 4*$y2 + $y3);
$s += $ds;
last if $ds < $precision;
$n -= 2;
}
return $s;
}