norm-quad-cdf.md

Running with SciPy:

#!/usr/bin/python3
import numpy as np
from scipy import integrate, stats

def norm_pdf(x):
    constant = 1.0 / np.sqrt(2*np.pi)
    return(constant * np.exp((-x**2) / 2.0) )

xs = [-1.0, 1.0]

# `quad()` takes scalar limits, so loop explicitly
res = [integrate.quad(norm_pdf, -np.inf, b, epsabs=1e-15, epsrel=0, limit=1000)[0]
       for b in xs]

print(res)
print(stats.norm.cdf(xs).tolist())  # vectorised
print(all(abs(r - c) < 1e-15 for r, c in zip(res, stats.norm.cdf(xs))))

Output:

[0.15865525393145707, 0.8413447460685431]
[0.15865525393145707, 0.8413447460685429]
True

---

Running in R:

#!/usr/bin/env Rscript -
xs <- c(-1, 1)

# `integrate()` takes scalar limits, so loop explicitly
res <- sapply(xs, function(b)
  integrate(dnorm, -Inf, b, abs.tol = 1e-15,
            rel.tol = 1e2 * .Machine$double.eps, subdivisions = 1000L)$value)

print(res, digits = 17)
print(pnorm(xs), digits = 17)  # vectorised
print(all(abs(res - pnorm(xs)) < 1e-15))

Output:

[1] 0.15865525393145707 0.84134474606854304
[1] 0.15865525393145705 0.84134474606854293
[1] TRUE

---

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];
# qagil: vectorised!
my ($res) = gslinteg_qagil( \&norm_pdf, $xvals , 1e-15, 0, 1000);

# CDF: vectorised
my $direct = gsl_cdf_gaussian_P( $xvals, 1);

say $res->string('%.18f');
say $direct->string('%.18f');
say all(abs( $res - $direct ) < 1e-15);

Output:

[0.158655253931457074 0.841344746068543148]
[0.158655253931457046 0.841344746068542926]
1