picomath

Scala

NormalCDFInverse.scala

import scala.math

//See http://www.johndcook.com/normal_cdf_inverse.html
class NormalCDFInverse {
//This is good for 0 < p \leq 0.5
//This is the rational approximation for the
//complementary cumulative distribution function
def rationalApproximation(t: Double) = {
// Abramowitz and Stegun formula 26.2.23.
// The absolute value of the error should be less than 4.5 e-4.
val c: Array[Double] = Array(2.515517, 0.802853, 0.010328)
val d: Array[Double] = Array(1.432788, 0.189269, 0.001308)
val numerator: Double = (c(2)*t + c(1))*t + c(0);
val denominator: Double = ((d(2)*t + d(1))*t + d(0))*t + 1.0;

t - numerator / denominator
}

//Only defined for 0 < p < 1
def normalCDFInverse(p: Double): Double = {
require(p > 0.0 && p < 1)
// See article above for explanation of this section.
if (p < 0.5) {
// F^-1(p) = - G^-1(p)
return -rationalApproximation( math.sqrt(-2.0*math.log(p)) );
} else {
// F^-1(p) = G^-1(1-p)
return rationalApproximation( math.sqrt(-2.0*math.log(1.0-p)) );
}
}
}

object NormalCDFInverse extends NormalCDFInverse