jetblack-options 0.1.0

use core::f64;

use libm::{fabs, log, sqrt};

fn r8poly_value(n: usize, a: &[f64], x: f64) -> f64 {
    let mut value = 0.0;

    for i in (0..n).rev() {
        value = value * x + a[i];
    }

    return value;
}

/// A function which computes the inverse of the Normal Cumulative Density
/// Function (CDF), by Michael Wichura;
///
/// This is a version of Applied Statistics Algorithm 241.
///
/// Original FORTRAN version:
///
/// Michael Wichura,
/// The Percentage Points of the Normal Distribution,
/// Algorithm AS 241,
/// Applied Statistics,
/// Volume 37, Number 3, pages 477-484, 1988.
///
fn r8_normal_01_cdf_inverse(p: f64) -> f64 {
    const A: [f64; 8] = [
        3.3871328727963666080,
        1.3314166789178437745e+2,
        1.9715909503065514427e+3,
        1.3731693765509461125e+4,
        4.5921953931549871457e+4,
        6.7265770927008700853e+4,
        3.3430575583588128105e+4,
        2.5090809287301226727e+3,
    ];
    const B: [f64; 8] = [
        1.0,
        4.2313330701600911252e+1,
        6.8718700749205790830e+2,
        5.3941960214247511077e+3,
        2.1213794301586595867e+4,
        3.9307895800092710610e+4,
        2.8729085735721942674e+4,
        5.2264952788528545610e+3,
    ];
    const C: [f64; 8] = [
        1.42343711074968357734,
        4.63033784615654529590,
        5.76949722146069140550,
        3.64784832476320460504,
        1.27045825245236838258,
        2.41780725177450611770e-1,
        2.27238449892691845833e-2,
        7.74545014278341407640e-4,
    ];
    const CONST1: f64 = 0.180625;
    const CONST2: f64 = 1.6;
    const D: [f64; 8] = [
        1.0,
        2.05319162663775882187,
        1.67638483018380384940,
        6.89767334985100004550e-1,
        1.48103976427480074590e-1,
        1.51986665636164571966e-2,
        5.47593808499534494600e-4,
        1.05075007164441684324e-9,
    ];
    const E: [f64; 8] = [
        6.65790464350110377720,
        5.46378491116411436990,
        1.78482653991729133580,
        2.96560571828504891230e-1,
        2.65321895265761230930e-2,
        1.24266094738807843860e-3,
        2.71155556874348757815e-5,
        2.01033439929228813265e-7,
    ];
    const F: [f64; 8] = [
        1.0,
        5.99832206555887937690e-1,
        1.36929880922735805310e-1,
        1.48753612908506148525e-2,
        7.86869131145613259100e-4,
        1.84631831751005468180e-5,
        1.42151175831644588870e-7,
        2.04426310338993978564e-15,
    ];
    const SPLIT1: f64 = 0.425;
    const SPLIT2: f64 = 5.0;

    if p <= 0.0 {
        return -1.0E+30;
    } else if 1.0 <= p {
        return 1.0E+30;
    }

    let q = p - 0.5;

    if fabs(q) <= SPLIT1 {
        let r = CONST1 - q * q;
        return q * r8poly_value(8, &A, r) / r8poly_value(8, &B, r);
    } else {
        let mut r = if q < 0.0 { p } else { 1.0 - p };

        if r <= 0.0 {
            return f64::NAN;
        }

        r = sqrt(-log(r));

        let mut value = if r <= SPLIT2 {
            r = r - CONST2;
            r8poly_value(8, &C, r) / r8poly_value(8, &D, r)
        } else {
            r = r - SPLIT2;
            r8poly_value(8, &E, r) / r8poly_value(8, &F, r)
        };

        if q < 0.0 {
            value = -value;
        }

        value
    }
}

pub fn inv_cdf(x: f64) -> f64 {
    r8_normal_01_cdf_inverse(x)
}

#[cfg(test)]
mod tests {
    use core::f64;

    use libm::fabs;

    use super::super::cdf::cdf;
    use super::*;

    fn is_close_to(actual: f64, expected: f64, threshold: f64) -> bool {
        let diff = fabs(actual - expected);
        diff < threshold
    }

    #[test]
    fn it_should_inverse_cdf() {
        for (x, fx, expected_x, expected_fx) in [
            (
                0.0000000000000000,
                0.5000000000000000,
                0.0000000000000000,
                0.5000000000000000,
            ),
            (
                0.1000000000000000,
                0.5398278372770290,
                0.10000000000000001,
                0.5398278372770290,
            ),
            (
                0.2000000000000000,
                0.5792597094391030,
                0.1999999999999999,
                0.5792597094391030,
            ),
            (
                0.3000000000000000,
                0.6179114221889526,
                0.29999999999999977,
                0.6179114221889526,
            ),
            (
                0.4000000000000000,
                0.6554217416103242,
                0.39999999999999997,
                0.6554217416103242,
            ),
            (
                0.5000000000000000,
                0.6914624612740131,
                0.49999999999999983,
                0.6914624612740131,
            ),
            (
                0.6000000000000000,
                0.7257468822499270,
                0.60000000000000164,
                0.72574688224992645,
            ),
            (
                0.7000000000000000,
                0.7580363477769270,
                0.69999999999999984,
                0.7580363477769270,
            ),
            (
                0.8000000000000000,
                0.7881446014166033,
                0.79999999999999982,
                0.7881446014166033,
            ),
            (
                0.9000000000000000,
                0.8159398746532405,
                0.9000000000000000,
                0.8159398746532405,
            ),
            (
                1.0000000000000000,
                0.8413447460685429,
                1.0000000000000000,
                0.8413447460685429,
            ),
            (
                1.5000000000000000,
                0.9331927987311419,
                1.4999999999999998,
                0.9331927987311419,
            ),
            (
                2.0000000000000000,
                0.9772498680518208,
                2.0000000000000000,
                0.9772498680518208,
            ),
            (
                2.5000000000000000,
                0.9937903346742240,
                2.5000000000000040,
                0.9937903346742240,
            ),
            (
                3.0000000000000000,
                0.9986501019683699,
                2.9999999999999969,
                0.9986501019683699,
            ),
            (
                3.5000000000000000,
                0.9997673709209645,
                3.4999999999999831,
                0.9997673709209645,
            ),
            (
                4.0000000000000000,
                0.9999683287581669,
                3.9999999999999996,
                0.9999683287581669,
            ),
        ] {
            let fx0 = cdf(x, 0.0, 1.0);
            assert!(is_close_to(fx0, expected_fx, f64::EPSILON));
            let x0 = inv_cdf(fx);
            assert!(is_close_to(x0, expected_x, f64::EPSILON));
        }
    }
}