com.quantarray.anaheim.numerics.StandardNormalDistribution.scala Maven / Gradle / Ivy
package com.quantarray.anaheim.numerics
/**
* Standard normal distribution.
*/
object StandardNormalDistribution
{
def pdf(x: Double): Double = NormalDistribution1.pdf(x, 0, 1)
/**
* Hart approximation of CDF.
*/
def cdf(x: Double): Double =
{
var n = 0.0
val absx = math.abs(x)
if (absx > 37)
{
n = 0
}
else
{
val exponential = math.exp(-absx * absx / 2)
if (absx < 7.07106781186547)
{
var build = 3.52624965998911E-02 * absx + 0.700383064443688
build = build * absx + 6.37396220353165
build = build * absx + 33.912866078383
build = build * absx + 112.079291497871
build = build * absx + 221.213596169931
build = build * absx + 220.206867912376
n = exponential * build
build = 8.83883476483184E-02 * absx + 1.75566716318264
build = build * absx + 16.064177579207
build = build * absx + 86.7807322029461
build = build * absx + 296.564248779674
build = build * absx + 637.333633378831
build = build * absx + 793.826512519948
build = build * absx + 440.413735824752
n = n / build
}
else
{
var build = absx + 0.65
build = absx + 4 / build
build = absx + 3 / build
build = absx + 2 / build
build = absx + 1 / build
n = exponential / build / 2.506628274631
}
}
if (x > 0)
n = 1 - n
n
}
/**
* Inverse of the cumulative distribution function.
* Equivalent of NORMSINV function in Excel.
*
* p: probability
*/
def cdfInverse(p: Double): Double =
{
if (p <= 0 || p >= 1)
{
throw new Exception("Number should be between 0 and 1")
}
/**
* Coefficients in rational approximations
*/
val a = List(-3.969683028665376e+01, 2.209460984245205e+02,
-2.759285104469687e+02, 1.383577518672690e+02,
-3.066479806614716e+01, 2.506628277459239e+00
)
val b = List(-5.447609879822406e+01, 1.615858368580409e+02,
-1.556989798598866e+02, 6.680131188771972e+01,
-1.328068155288572e+01
)
val c = List(-7.784894002430293e-03, -3.223964580411365e-01,
-2.400758277161838e+00, -2.549732539343734e+00,
4.374664141464968e+00, 2.938163982698783e+00
)
val d = List(7.784695709041462e-03, 3.224671290700398e-01,
2.445134137142996e+00, 3.754408661907416e+00
)
/**
* Define break-points.
*/
val plow = 0.02425
val phigh = 1 - plow
/**
* Rational approximation for lower region:
*/
if (p < plow)
{
val q = math.sqrt(-2 * math.log(p))
return (((((c.head * q + c(1)) * q + c(2)) * q + c(3)) * q + c(4)) * q + c(5)) /
((((d.head * q + d(1)) * q + d(2)) * q + d(3)) * q + 1)
}
/**
* Rational approximation for upper region:
*/
if (phigh < p)
{
val q = math.sqrt(-2 * math.log(1 - p))
return -(((((c.head * q + c(1)) * q + c(2)) * q + c(3)) * q + c(4)) * q + c(5)) /
((((d.head * q + d(1)) * q + d(2)) * q + d(3)) * q + 1)
}
/**
* Rational approximation for central region:
*/
val q = p - 0.5
val r = q * q
(((((a.head * r + a(1)) * r + a(2)) * r + a(3)) * r + a(4)) * r + a(5)) * q /
(((((b.head * r + b(1)) * r + b(2)) * r + b(3)) * r + b(4)) * r + 1)
}
}