// Hypothesis Testing
// Martin Sewell
// 30 May 2020
// Functions gathered from the web, tested and selected for replicating Excel functions to the highest degree of accuracy.

//Equivalent Excel functions
//NORM.DIST(x, mean, standard_dev, TRUE)
//NORM.DIST(x, mean, standard_dev, FALSE)
//NORM.INV(probability, mean, standard_dev)
//NORM.S.DIST(z, TRUE)
//NORM.S.DIST(z, FALSE)
//NORM.S.INV(probability)

#define _USE_MATH_DEFINES

#include <iostream>


//NORM.DIST(x, mean, standard_dev, TRUE)
const long double eps = std::numeric_limits<long double>::epsilon();
long double NormDistCDF(long double x, long double mu, long double sigma)
{
    // This algorithm is ported from dcdflib:
    // Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
    // Package of Special Function Routines and Test Drivers"
    // acm Transactions on Mathematical Software. 19, 22-32.
    
    int i;
    long double del, xden, xnum, xsq;
    long double result, ccum;
    long double arg = (x - mu) / sigma;
    const long double sixten = 1.60e0;
    const long double sqrpi = 3.9894228040143267794e-1;
    const long double thrsh = 0.66291e0;
    const long double root32 = 5.656854248e0;
    const long double zero = 0.0e0;
    const long double min = eps; // Double.Epsilon;
    long double z = arg;
    long double y = abs(z);
    const long double half = 0.5e0;
    const long double one = 1.0e0;

    if (sigma < 0)
    {
        throw "The standard deviation must be positive.";
    }
    long double a[] =
    {
        2.2352520354606839287e00, 1.6102823106855587881e02, 1.0676894854603709582e03,
        1.8154981253343561249e04, 6.5682337918207449113e-2
    };
    long double b[] =
    {
        4.7202581904688241870e01, 9.7609855173777669322e02, 1.0260932208618978205e04,
        4.5507789335026729956e04
    };
    long double c[] =
    {
        3.9894151208813466764e-1, 8.8831497943883759412e00, 9.3506656132177855979e01,
        5.9727027639480026226e02, 2.4945375852903726711e03, 6.8481904505362823326e03,
        1.1602651437647350124e04, 9.8427148383839780218e03, 1.0765576773720192317e-8
    };
    long double d[] =
    {
        2.2266688044328115691e01, 2.3538790178262499861e02, 1.5193775994075548050e03,
        6.4855582982667607550e03, 1.8615571640885098091e04, 3.4900952721145977266e04,
        3.8912003286093271411e04, 1.9685429676859990727e04
    };
    long double p[] =
    {
        2.1589853405795699e-1, 1.274011611602473639e-1, 2.2235277870649807e-2,
        1.421619193227893466e-3, 2.9112874951168792e-5, 2.307344176494017303e-2
    };
    long double q[] =
    {
        1.28426009614491121e00, 4.68238212480865118e-1, 6.59881378689285515e-2,
        3.78239633202758244e-3, 7.29751555083966205e-5
    };
    if (y <= thrsh)
    {
        //
        // Evaluate  anorm  for  |X| <= 0.66291
        //
        xsq = zero;
        if (y > std::numeric_limits<long double>::epsilon()) xsq = z * z;
        xnum = a[4] * xsq;
        xden = xsq;
        for (i = 0; i < 3; i++) {
            xnum = (xnum + a[i]) * xsq;
            xden = (xden + b[i]) * xsq;
        }
        result = z * (xnum + a[3]) / (xden + b[3]);
        long double temp = result;
        result = half + temp;
    }
    else if (y <= root32)
    {
        xnum = c[8] * y;
        xden = y;
        for (i = 0; i < 7; i++) {
            xnum = (xnum + c[i]) * y;
            xden = (xden + d[i]) * y;
        }
        result = (xnum + c[7]) / (xden + d[7]);
        xsq = floor(y * sixten) / sixten;
        del = (y - xsq) * (y + xsq);
        result = exp(-(xsq * xsq * half)) * exp(-(del * half)) * result;
        ccum = one - result;
        if (z > zero) {
            result = ccum;
        }
    }
    else
    {
        xsq = one / (z * z);
        xnum = p[5] * xsq;
        xden = xsq;
        for (i = 0; i < 4; i++)
        {
            xnum = (xnum + p[i]) * xsq;
            xden = (xden + q[i]) * xsq;
        }
        result = xsq * (xnum + p[4]) / (xden + q[4]);
        result = (sqrpi - result) / y;
        xsq = floor(z * sixten) / sixten;
        del = (z - xsq) * (z + xsq);
        result = exp(-(xsq * xsq * half)) * exp(-(del * half)) * result;
        ccum = one - result;
        if (z > zero) {
            result = ccum;
        }
    }

    if (result < min)
        result = 0.0e0;
    return result;
}


//NORM.DIST(x, mean, standard_dev, FALSE)
long double NormDistPDF(long double x, long double mu, long double sigma)
{
    if (sigma < 0)
    {
        throw "The standard deviation must be positive.";
    }
    long double en = -(x - mu) * (x - mu);
    long double ed = 2 * sigma * sigma;
    long double e = en / ed;
    long double d = sigma * sqrt(2 * M_PI);

    return (1 / d) * exp(e);
}


//NORM.INV(probability, mean, standard_dev)
long double NormInv(long double p, long double mu, long double sigma)
{
    if (p < 0 || p > 1)
    {
        throw "The probability must be greater than 0 and less than 1.";
    }
    if (sigma < 0)
    {
        throw "The standard deviation must be positive.";
    }
    if (p == 0)
    {
        return -999999; //-Infinity;
    }
    if (p == 1)
    {
        return 999999; //Infinity;
    }
    if (sigma == 0)
    {
        return mu;
    }

    long double q, r, val;
    q = p - 0.5;

    /*-- use AS 241 --- */
    /* double ppnd16_(double *p, long *ifault)*/
    /*      ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3
            Produces the normal deviate Z corresponding to a given lower
            tail area of P; Z is accurate to about 1 part in 10**16.
    */
    if (abs(q) <= .425)
    {/* 0.075 <= p <= 0.925 */
        r = .180625 - q * q;
        val =
            q * (((((((r * 2509.0809287301226727 +
                33430.575583588128105) * r + 67265.770927008700853) * r +
                45921.953931549871457) * r + 13731.693765509461125) * r +
                1971.5909503065514427) * r + 133.14166789178437745) * r +
                3.387132872796366608)
            / (((((((r * 5226.495278852854561 +
                28729.085735721942674) * r + 39307.89580009271061) * r +
                21213.794301586595867) * r + 5394.1960214247511077) * r +
                687.1870074920579083) * r + 42.313330701600911252) * r + 1);
    }
    else
    { /* closer than 0.075 from {0,1} boundary */
        /* r = min(p, 1-p) < 0.075 */
        if (q > 0)
            r = 1 - p;
        else
            r = p;
        r = sqrt(-log(r));
        /* r = sqrt(-log(r))  <==>  min(p, 1-p) = exp( - r^2 ) */
        if (r <= 5)
        { /* <==> min(p,1-p) >= exp(-25) ~= 1.3888e-11 */
            r += -1.6;
            val = (((((((r * 7.7454501427834140764e-4 +
                .0227238449892691845833) * r + .24178072517745061177) *
                r + 1.27045825245236838258) * r +
                3.64784832476320460504) * r + 5.7694972214606914055) *
                r + 4.6303378461565452959) * r +
                1.42343711074968357734)
                / (((((((r *
                    1.05075007164441684324e-9 + 5.475938084995344946e-4) *
                    r + .0151986665636164571966) * r +
                    .14810397642748007459) * r + .68976733498510000455) *
                    r + 1.6763848301838038494) * r +
                    2.05319162663775882187) * r + 1);
        }
        else
        { /* very close to  0 or 1 */
            r += -5;
            val = (((((((r * 2.01033439929228813265e-7 +
                2.71155556874348757815e-5) * r +
                .0012426609473880784386) * r + .026532189526576123093) *
                r + .29656057182850489123) * r +
                1.7848265399172913358) * r + 5.4637849111641143699) *
                r + 6.6579046435011037772)
                / (((((((r *
                    2.04426310338993978564e-15 + 1.4215117583164458887e-7) *
                    r + 1.8463183175100546818e-5) * r +
                    7.868691311456132591e-4) * r + .0148753612908506148525)
                    * r + .13692988092273580531) * r +
                    .59983220655588793769) * r + 1);
        }
        if (q < 0.0)
        {
            val = -val;
        }
    }
    return mu + sigma * val;
}


//NORM.S.DIST(z, TRUE)
long double NormSDistCDF(long double z)
{
    return 0.5 * erfc(-z * M_SQRT1_2);
}


//NORM.S.DIST(z, FALSE)
long double NormSDistPDF(long double z)
{
    long double en = -z * z;
    long double ed = 2;
    long double e = en / ed;
    long double d = sqrt(2 * M_PI);

    return (1 / d) * exp(e);
}


//NORM.S.INV(probability)
#define  A1  (-3.969683028665376e+01)
#define  A2   2.209460984245205e+02
#define  A3  (-2.759285104469687e+02)
#define  A4   1.383577518672690e+02
#define  A5  (-3.066479806614716e+01)
#define  A6   2.506628277459239e+00

#define  B1  (-5.447609879822406e+01)
#define  B2   1.615858368580409e+02
#define  B3  (-1.556989798598866e+02)
#define  B4   6.680131188771972e+01
#define  B5  (-1.328068155288572e+01)

#define  C1  (-7.784894002430293e-03)
#define  C2  (-3.223964580411365e-01)
#define  C3  (-2.400758277161838e+00)
#define  C4  (-2.549732539343734e+00)
#define  C5   4.374664141464968e+00
#define  C6   2.938163982698783e+00

#define  D1   7.784695709041462e-03
#define  D2   3.224671290700398e-01
#define  D3   2.445134137142996e+00
#define  D4   3.754408661907416e+00

#define P_LOW   0.02425
/* P_high = 1 - p_low*/
#define P_HIGH  0.97575

long double NormSInv(long double p)
{
    if (p < 0 || p > 1)
    {
        throw "The probability must be greater than 0 and less than 1.";
    }
    long double x = 0.0;
    long double q, r, u, e;
    if ((0 < p) && (p < P_LOW)) {
        q = sqrt(-2 * log(p));
        x = (((((C1 * q + C2) * q + C3) * q + C4) * q + C5) * q + C6) / ((((D1 * q + D2) * q + D3) * q + D4) * q + 1);
    }
    else {
        if ((P_LOW <= p) && (p <= P_HIGH)) {
            q = p - 0.5;
            r = q * q;
            x = (((((A1 * r + A2) * r + A3) * r + A4) * r + A5) * r + A6) * q / (((((B1 * r + B2) * r + B3) * r + B4) * r + B5) * r + 1);
        }
        else {
            if ((P_HIGH < p) && (p < 1)) {
                q = sqrt(-2 * log(1 - p));
                x = -(((((C1 * q + C2) * q + C3) * q + C4) * q + C5) * q + C6) / ((((D1 * q + D2) * q + D3) * q + D4) * q + 1);
            }
        }
    }
    if ((0 < p) && (p < 1)) {
        e = 0.5 * erfc(-x / sqrt(2)) - p;
        u = e * sqrt(2 * M_PI) * exp(x * x / 2);
        x = x - u / (1 + x * u / 2);
    }

    return x;
}



//Equivalent Excel functions
//NORM.DIST(x, mean, standard_dev, TRUE)
//NORM.DIST(x, mean, standard_dev, FALSE)
//NORM.INV(probability, mean, standard_dev)
//NORM.S.DIST(z, TRUE)
//NORM.S.DIST(z, FALSE)
//NORM.S.INV(probability)

int main()
{
    std::cout << NormDistCDF(1, 0, 1) << std::endl; // NORM.DIST(1, 0, 1, TRUE)
    std::cout << NormDistPDF(1, 0, 1) << std::endl; // NORM.DIST(1, 0, 1, FALSE)
    std::cout << NormInv(0.3, 0, 1) << std::endl;   // NORM.INV(0.3, 0, 1)
    std::cout << NormSDistCDF(1) << std::endl;      // NORM.S.DIST(1, TRUE)
    std::cout << NormSDistPDF(1) << std::endl;      // NORM.S.DIST(1, FALSE)
    std::cout << NormSInv(0.3) << std::endl;        // NORM.S.INV(0.3)

    std::cin.get();
    return 0;
}