Skip to content

Instantly share code, notes, and snippets.

@gusty
Last active November 17, 2022 05:15
Show Gist options
  • Save gusty/06ab1eeaba3ce80dcfcc to your computer and use it in GitHub Desktop.
Save gusty/06ab1eeaba3ce80dcfcc to your computer and use it in GitHub Desktop.
Generic CDF
#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
*)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment