normal-dist.dm42

// Normal distribution functions
// 
// Im trying to only use my dm42 for my statistics class (bad idea probally),
// hoping we don't end up needing any more complacated functions...

// Refrences:
// -  Winitzki, Sergei. "A handy approximation for the error function and its inverse." A lecture note obtained through private communication (2008).
// - https://en.wikipedia.org/wiki/Cumulative_distribution_function
// - https://stackoverflow.com/questions/457408/is-there-an-easily-available-implementation-of-erf-for-python

// Normal cumulative distribution function
//
// (1 + erf((x - mu) / sigma / sqrt(2))) / 2
export def cdf {
    MVAR "X"
    MVAR "MU"
    MVAR "SIGMA"
    VARMENU cdf
    STOP

    RCL "X"
    RCL "MU"
    -
    RCL "SIGMA"
    /
    2
    SQRT
    /
    erf()
    1, +
    2, /

    exit()
}

// Inverse normal cumulative distribution function
//
// mu + sigma * sqrt(2) * inv_erf(2 * x - 1)
export def icdf {
    MVAR "X"
    MVAR "MU"
    MVAR "SIGMA"
    VARMENU icdf
    STOP

    2
    RCL "X"
    *
    1
    -
    inv_erf()
    2
    SQRT
    *
    RCL "SIGMA"
    *
    RCL "MU"
    +

    exit()
}

def exit {
    CLMENU
    if { FS? 03 } {
        SF 27
    }
}

// Error function
//
// sign = signum(x)
// x = abs(x)
// 
// a1 =  0.254829592
// a2 = -0.284496736
// a3 =  1.421413741
// a4 = -1.453152027
// a5 =  1.061405429
// p  =  0.3275911
// 
// t = 1.0/(1.0 + p*x)
// y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1) * t * e^(-x^2)
// sign*y
def erf {
    // Take the absolute value of x
    // also save the sign of x for later
    DUP
    ABS
    STO 00 // |X| in 00
    DROP
    SIGN
    STO 01 // sign(X) in 01
    DROP

    // Calculate t
    0.3275911 // p
    RCL 00
    *
    1,+
    1/X
    STO 02 // t in 02
    DROP

    // Calculate y
    1
    1.061405429 // a5
    RCL 02, *
    -1.453152027 // a4
    +
    RCL 02, *
    1.421413741 // a3
    +
    RCL 02, *
    -0.284496736 // a2
    +
    RCL 02, *
    0.254829592 // a1
    +
    RCL 02, *

    // * e^(-x^2)
    RCL 00
    X↑2
    +/-
    E↑X
    *
    -

    // Return sign
    RCL 01
    *
}

// Inverse error function
//
// s = signum(x)
// x = abs(x)
//
// a = (8 * (pi - 3)) / (3 * pi * (4 - pi)) # ~0.140012
// s * sqrt(sqrt((2 / (pi * a) + (log(1 - x ** 2)) / 2) ** 2 - log(1 - x ** 2) / a) - (2 / (pi * a) + log(1 - x ** 2) / 2))
def inv_erf {
    DUP
    ABS
    STO 00 // |X| in 00
    DROP
    SIGN
    STO 01 // sign(X) in 01
    DROP

    // Calculate of a
    8
    PI
    3
    -
    *
    4
    PI
    -
    PI
    *
    3
    *
    /
    STO 02 // a in 02
    DROP

    // Calculate
    2
    PI
    RCL 02
    *
    /
    1
    RCL 00
    X↑2
    -
    LN
    2
    /
    +
    X↑2
    1
    RCL 00
    X↑2
    -
    LN
    RCL 02
    /
    -
    SQRT
    2
    PI
    RCL 02
    *
    /
    1
    RCL 00
    X↑2
    -
    LN
    2
    /
    +
    -
    SQRT

    // Restore sign
    RCL 01
    *
}