! Copyright (C) 2011 John Benediktsson ! See https://factorcode.org/license.txt for BSD license USING: combinators kernel locals math math.constants math.functions sequences ; IN: picomath ! Standalone error function erf(x) ! https://www.johndcook.com/blog/2009/01/19/stand-alone-error-function-erf/ :: erf ( x -- value ) x 0 >= 1 -1 ? :> sign x abs :> x! p x * 1 + recip :> t a5 t * a4 + t * a3 + t * a2 + t * a1 + t * x x neg * e^ * 1 swap - :> y sign y * ; :: expm1 ( x -- value ) x abs 1e-5 < [ x sq 0.5 * x + ] [ x e^ 1.0 - ] if ; ! Standalone implementation of phi(x) :: phi ( x -- value ) x 0 >= 1 -1 ? :> sign x abs 2 sqrt / :> x! p x * 1 + recip :> t a5 t * a4 + t * a3 + t * a2 + t * a1 + t * x x neg * e^ * 1 swap - :> y sign y * 1 + 2 / ; : log-factorial ( n -- value ) { { [ dup 0 < ] [ "invalid input" throw ] } { [ dup 254 > ] [ 1 + [| x | x 0.5 - x log * x - 0.5 2 pi * log * + 1.0 12.0 x * / + ] call ] } [ lf nth ] } cond ; :: log-one-plus-x ( x -- value ) x -1.0 <= [ "argument must be > -1" throw ] when x abs 1e-4 > [ 1.0 x + log ] [ -0.5 x * 1.0 + x * ] if ; :: rational-approximation ( t -- value ) c2 t * c1 + t * c0 + :> numerator d2 t * d1 + t * d0 + t * 1.0 + :> denominator t numerator denominator / - ; :: normal-cdf-inverse ( p -- value ) p [ 0 > ] [ 1 < ] bi and [ p throw ] unless p 0.5 < [ p log -2.0 * sqrt rational-approximation neg ] [ p 1.0 - log -2.0 * sqrt rational-approximation ] if ; DEFER: gamma :: log-gamma ( x -- value ) x 0 <= [ "Invalid input" throw ] when x 12 < [ x gamma abs log ] [ 1.0 x x * / :> z 7 c nth 7 reverse [ [ z * ] [ c nth ] bi* + ] each x / :> series x 0.5 - x log * x - halfLogTwoPi + series + ] if ; y! 0 :> n! y 1.0 < :> arg-was-less-than-one arg-was-less-than-one [ y 1.0 + y! ] [ y floor >integer 1 - n! y n - y! ] if 0.0 :> num! 1.0 :> den! y 1 - :> z! 8 [ [ P nth num + z * num! ] [ Q nth den z * + den! ] bi ] each num den / 1.0 + arg-was-less-than-one [ y 1.0 - / ] [ n [ y * y 1.0 + y! ] times ] if ; PRIVATE> :: gamma ( x -- value ) x 0 <= [ "Invalid input" throw ] when x { { [ dup 0.001 < ] [ GAMMA * 1.0 + x * 1.0 swap / ] } { [ dup 12.0 < ] [ (gamma) ] } { [ dup 171.624 > ] [ drop 1/0. ] } [ log-gamma e^ ] } cond ;