math.extras: adding sterling numbers, bernoulli numbers, chi-square and inverse-chi-square words.

extra/math/extras/extras-docs.factor

@@ -0,0 +1,23 @@
+! Copyright (C) 2012 John Benediktsson
+! See http://factorcode.org/license.txt for BSD license
+
+USING: help.markup help.syntax math sequences ;
+
+IN: math.extras
+
+HELP: bernoulli
+{ $values { "p" integer } { "n" rational } }
+{ $description "Return the Bernoulli number " { $snippet "p" } "." } ;
+
+HELP: sterling
+{ $values { "n" integer } { "k" integer } { "x" integer } }
+{ $description "Return the Stirling number of the second kind for a set with " { $snippet "n" } " elements partitioned into " { $snippet "k" } " disjoint non-empty sets." } ;
+
+HELP: chi2
+{ $values { "actual" sequence } { "expected" sequence } { "n" real } }
+{ $description "Return the chi-squared metric between " { $snippet "actual" } " and " { $snippet "expected" } " observations." } ;
+
+HELP: chi2P
+{ $values { "chi" real } { "df" real } { "p" real } }
+{ $description "Returns the inverse chi-squared value according to " { $snippet "P(chi|df) = P(df/2,chi/2)" } "." } ;
+

extra/math/extras/extras.factor

@@ -0,0 +1,53 @@
+! Copyright (C) 2012 John Benediktsson
+! See http://factorcode.org/license.txt for BSD license
+
+USING: combinators.short-circuit kernel math math.combinatorics
+math.functions math.order math.ranges math.statistics
+math.vectors memoize sequences ;
+
+IN: math.extras
+
+<PRIVATE
+
+DEFER: sterling
+
+: (sterling) ( n k -- x )
+    [ [ 1 - ] bi@ sterling ]
+    [ [ 1 - ] dip sterling ]
+    [ nip * + ] 2tri ;
+
+PRIVATE>
+
+MEMO: sterling ( n k -- x )
+    2dup { [ = ] [ nip 1 = ] } 2||
+    [ 2drop 1 ] [ (sterling) ] if ;
+
+<PRIVATE
+
+DEFER: bernoulli
+
+: (bernoulli) ( p -- n )
+    [ iota ] [ 1 + ] bi [
+        0 [ [ nCk ] [ bernoulli * ] bi + ] with reduce
+    ] keep recip neg * ;
+
+PRIVATE>
+
+MEMO: bernoulli ( p -- n )
+    [ 1 ] [ (bernoulli) ] if-zero ;
+
+: chi2 ( actual expected -- n )
+    0 [ dup 0 > [ [ - sq ] keep / + ] [ 2drop ] if ] 2reduce ;
+
+<PRIVATE
+
+: df-check ( df -- )
+    even? [ "odd degrees of freedom" throw ] unless ;
+
+: (chi2P) ( chi/2 df/2 -- p )
+    [1,b) dupd n/v cum-product swap neg exp [ v*n sum ] keep + ;
+
+PRIVATE>
+
+: chi2P ( chi df -- p )
+    dup df-check [ 2.0 / ] [ 2 /i ] bi* (chi2P) 1.0 min ;