jetblack-options 0.2.1

//! The cumulative bivariate normal distribution function

use core::f64::consts::PI;

use libm::{asin, exp, fabs, fmax, sin, sqrt};

use super::cnd::CND;

fn sqr(x: f64) -> f64 {
    x * x
}

#[allow(non_snake_case)]
struct Factor {
    W: &'static [f64],
    XX: &'static [f64],
}
const F3: Factor = Factor {
    W: &[0.17132449237917, 0.360761573048138, 0.46791393457269],
    XX: &[-0.932469514203152, -0.661209386466265, -0.238619186083197],
};
const F6: Factor = Factor {
    W: &[
        4.71753363865118E-02,
        0.106939325995318,
        0.160078328543346,
        0.203167426723066,
        0.233492536538355,
        0.249147045813403,
    ],
    XX: &[
        -0.981560634246719,
        -0.904117256370475,
        -0.769902674194305,
        -0.587317954286617,
        -0.36783149899818,
        -0.125233408511469,
    ],
};
const F10: Factor = Factor {
    W: &[
        1.76140071391521E-02,
        4.06014298003869E-02,
        6.26720483341091E-02,
        8.32767415767048E-02,
        0.10193011981724,
        0.118194531961518,
        0.131688638449177,
        0.142096109318382,
        0.149172986472604,
        0.152753387130726,
    ],
    XX: &[
        -0.993128599185095,
        -0.963971927277914,
        -0.912234428251326,
        -0.839116971822219,
        -0.746331906460151,
        -0.636053680726515,
        -0.510867001950827,
        -0.37370608871542,
        -0.227785851141645,
        -7.65265211334973E-02,
    ],
};

///   A function for computing bivariate normal probabilities.
///     Alan Genz
///     Department of Mathematics
///     Washington State University
///     Pullman, WA 99164-3113
///     Email : alangenz@wsu.edu
///  This function is based on the method described by
///      Drezner, Z and G.O. Wesolowsky, (1990),
///      On the computation of the bivariate normal integral,
///      Journal of Statist. Comput. Simul. 35, pp. 101-107,
///  with major modifications for double precision, and for |R| close to 1.
/// This code was originally translated into VBA by Graeme West
#[allow(non_snake_case)]
pub fn cbnd(x: f64, y: f64, rho: f64) -> f64 {
    let Factor { W, XX } = if fabs(rho) < 0.3 {
        &F3
    } else if fabs(rho) < 0.75 {
        &F6
    } else {
        &F10
    };

    let h = -x;
    let mut k = -y;
    let mut hk = h * k;
    let mut BVN = 0.0;

    if fabs(rho) < 0.925 {
        if fabs(rho) > 0.0 {
            let hs = (h * h + k * k) / 2.0;
            let asr = asin(rho);
            for i in 0..W.len() {
                for ISs in [-1.0, 1.0] {
                    let sn = sin(asr * (ISs * XX[i] + 1.0) / 2.0);
                    BVN = BVN + W[i] * exp((sn * hk - hs) / (1.0 - sn * sn));
                }
            }
            BVN = BVN * asr / (4.0 * PI);
        }
        BVN = BVN + CND(-h) * CND(-k)
    } else {
        if rho < 0.0 {
            k = -k;
            hk = -hk;
        }
        if fabs(rho) < 1.0 {
            let Ass = (1.0 - rho) * (1.0 + rho);
            let A = sqrt(Ass);
            let bs = (h - k) * (h - k);
            let c = (4.0 - hk) / 8.0;
            let d = (12.0 - hk) / 16.0;
            let asr = -(bs / Ass + hk) / 2.0;
            if asr > -100.0 {
                BVN = A
                    * exp(asr)
                    * (1.0 - c * (bs - Ass) * (1.0 - d * bs / 5.0) / 3.0 + c * d * Ass * Ass / 5.0);
            }
            if -hk < 100.0 {
                let b = sqrt(bs);
                BVN = BVN
                    - exp(-hk / 2.0)
                        * sqrt(2.0 * PI)
                        * CND(-b / A)
                        * b
                        * (1.0 - c * bs * (1.0 - d * bs / 5.0) / 3.0);
            }
            let A = A / 2.0;
            for i in 0..W.len() {
                for ISs in [-1.0, 1.0] {
                    let xs = sqr(A * (ISs * XX[i] + 1.0));
                    let rs = sqrt(1.0 - xs);
                    let asr = -(bs / xs + hk) / 2.0;
                    if asr > -100.0 {
                        BVN = BVN
                            + A * W[i]
                                * exp(asr)
                                * (exp(-hk * (1.0 - rs) / (2.0 * (1.0 + rs))) / rs
                                    - (1.0 + c * xs * (1.0 + d * xs)));
                    }
                }
            }
            BVN = -BVN / (2.0 * PI);
        }
        if rho > 0.0 {
            BVN = BVN + CND(-fmax(h, k));
        } else {
            BVN = -BVN;
            if k > h {
                BVN = BVN + CND(k) - CND(h);
            }
        }
    }

    BVN
}

#[cfg(test)]
mod tests {
    use libm::fabs;

    use super::*;

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

    #[test]
    fn it_should_calc_cbnd() {
        for (x, y, rho, expected, threshold) in [
            (-2.000, 1.000, 0.500, 0.02260327218216498, 1e-10),
            (-1.000, 1.000, 0.500, 0.1548729518586029, 1e-10),
            (0.000, 1.000, 0.500, 0.4687429526451681, 1e-10),
            (1.000, 1.000, 0.500, 0.7452035868467498, 1e-10),
            (2.000, 1.000, 0.500, 0.8318608311308804, 1e-10),
            (3.000, 1.000, 0.500, 0.8410314268855684, 1e-10),
            (-0.200, 0.500, -0.900, 0.1377019384919554, 1e-10),
            (-0.200, 0.500, -0.800, 0.1621749501739237, 1e-10),
            (-0.200, 0.500, -0.700, 0.182741124323837, 1e-10),
            (-0.200, 0.500, -0.600, 0.2010067421350173, 1e-10),
            (-0.200, 0.500, -0.500, 0.2177751156457083, 1e-10),
            (-0.200, 0.500, -0.400, 0.2335088436667788, 1e-10),
            (-0.200, 0.500, -0.300, 0.2485057781836367, 1e-10),
            (-0.200, 0.500, -0.200, 0.2629747825156674, 1e-10),
            (-0.200, 0.500, -0.100, 0.277072982340618, 1e-10),
            (-0.200, 0.500, 0.000, 0.2909261168683812, 1e-10),
            (-0.200, 0.500, 0.100, 0.3046406378728856, 1e-10),
            (-0.200, 0.500, 0.200, 0.318311355328648, 1e-10),
            (-0.200, 0.500, 0.300, 0.3320262544201822, 1e-10),
            (-0.200, 0.500, 0.400, 0.3458686754635103, 1e-10),
            (-0.200, 0.500, 0.500, 0.3599150462316804, 1e-10),
            (-0.200, 0.500, 0.600, 0.3742210899873292, 1e-10),
            (-0.200, 0.500, 0.700, 0.38877064052855, 1e-10),
            (-0.200, 0.500, 0.800, 0.4032765198363306, 1e-10),
            (-0.200, 0.500, 0.900, 0.4162100291976345, 1e-10),
            (1.000, 0.500, 0.673, 0.6508271498843442, 1e-10),
            (2.000, 1.000, 0.500, 0.8318608311308804, 1e-10),
            (0.000, 0.000, -1.000, 0.0, 1e-10),
            (0.000, 0.000, -0.500, 0.1666666666666667, 1e-10),
            (0.000, 0.000, 0.000, 0.25, 1e-10),
            (0.000, 0.000, 0.500, 0.3333333333333333, 1e-10),
            (0.000, 0.000, 1.000, 0.5, 1e-10),
            (1.000, 1.000, 0.500, 0.7452035868467498, 1e-10),
            (1.000, -1.000, 0.500, 0.1548729518586028, 1e-10),
            (-1.000, 1.000, 0.500, 0.1548729518586029, 1e-10),
            (-1.000, -1.000, 0.500, 0.06251409470966388, 1e-10),
            (1.000, 1.000, 0.500, 0.7452035868467498, 1e-10),
            (1.000, -1.000, 0.500, 0.1548729518586028, 1e-10),
            (-1.000, 1.000, 0.500, 0.1548729518586029, 1e-10),
            (-1.000, -1.000, 0.500, 0.06251409470966388, 1e-10),
            (0.707, 0.707, 0.500, 0.6337020457780796, 1e-4),
        ] {
            let actual = cbnd(x, y, rho);
            assert!(is_close_to(actual, expected, threshold));
        }
    }
}