Generic CDF

genericCDF.fsx

#r "nuget: FSharpPlus,1.2"

open FSharpPlus
open FSharpPlus.Math.Generic

type Ratio =
    struct
        val Numerator   : bigint
        val Denominator : bigint
        new (numerator: bigint, denominator: bigint) = {Numerator = numerator; Denominator = denominator}
    end
    override this.ToString() = this.Numerator.ToString() + " % " + this.Denominator.ToString()

let inline ratio (a: bigint) (b: bigint) : Ratio =
    if b = 0I then failwith "Ratio.%: zero denominator"
    let a, b = if b < 0I then (-a, -b) else (a, b)
    let gcd = gcd a b
    Ratio (a / gcd, b / gcd)

let Ratio (x,y) = x </ratio/> y

type Ratio with
    static member (/) (a:Ratio, b:Ratio) = (a.Numerator * b.Denominator) </ratio/> (a.Denominator * b.Numerator)
    static member (+) (a:Ratio, b:Ratio) = (a.Numerator * b.Denominator + b.Numerator * a.Denominator) </ratio/> (a.Denominator * b.Denominator)
    static member (-) (a:Ratio, b:Ratio) = (a.Numerator * b.Denominator - b.Numerator * a.Denominator) </ratio/> (a.Denominator * b.Denominator)
    static member (*) (a:Ratio, b:Ratio) = (a.Numerator * b.Numerator) </ratio/> (a.Denominator * b.Denominator)

    static member Abs        (r:Ratio) = (abs    r.Numerator) </ratio/> r.Denominator
    static member Signum     (r:Ratio) = (signum r.Numerator) </ratio/> 1I
    static member FromBigInt (x:bigint) = fromBigInt x </ratio/> 1I
    static member (~-)       (r:Ratio) = -(r.Numerator) </ratio/> r.Denominator

let inline fromRational (x:Ratio) : 'Fractional =
    let f x = explicit x : 'Fractional
    f ( x.Numerator) / f ( x.Denominator)



let inline CDF (x:^T) : ^T = 
    let num x = fromRational (x </ratio/> 1000000000I)
    let (b1,b2,b3)  = (num 319381530I   , num -356563782I  , num 1781477937I)
    let (b4,b5)     = (num -1821255978I , num 1330274429I)
    let (p , c )    = (num  0231641900I , num 0398942280I)
    let (zero, one, two) = 0G, 1G, 2G
    if x > zero then
        let t = one / (one + p * x) 
        (one - c * exp( -x * x / two)* t * (t*(t*(t*(t*b5+b4)+b3)+b2)+b1)) 
    else
        let t = one / (one - p * x) 
        (c * exp( -x * x / two)* t * (t*(t*(t*(t*b5+b4)+b3)+b2)+b1))

(*
> CDF (2.3) ;;
val it: float = 0.989275919

> CDF (2.3f) ;;
val it: float32 = 0.9892759323f
*)