Running with SciPy:
#!/usr/bin/python3
import scipy.integrate as integ
import numpy as np
def normalProbabilityDensity(x):
constant = 1.0 / np.sqrt(2*np.pi)
return(constant * np.exp((-x**2) / 2.0) )
neg_percentile, _ = integ.quad(normalProbabilityDensity, np.NINF, -1.00)
pos_percentile, _ = integ.quad(normalProbabilityDensity, np.NINF, 1.00)
print('Neg: ', neg_percentile)
print('Pos: ', pos_percentile)
# Neg: 0.15865525393145707
# Pos: 0.8413447460685435
Running in R:
#!/usr/bin/env Rscript -
pnorm(c(-1,1))
# [1] 0.1586553 0.8413447
Running with PDL:
#!/usr/bin/env perl
use PDL;
use PDL::Constants qw(PI);
use PDL::GSL::INTEG;
use PDL::GSL::CDF;
use feature qw(state say);
sub norm_pdf {
my ($x) = @_;
state $constant = 1 / sqrt(2*PI);
return $constant * exp(- $x**2 / 2);
}
my $xvals = [-1, 1];
my ($res) = gslinteg_qagil( \&norm_pdf, $xvals , 1e-15, 0, 1000);
my $direct = gsl_cdf_gaussian_P( $xvals, 1);
say $res->string('%.18f');
say $direct->string('%.18f');
say all(abs( $res - $direct ) < 1e-15);
# [0.158655253931457074 0.841344746068543148]
# [0.158655253931457046 0.841344746068542926]
# 1