1 ! Copyright (C) 2012 John Benediktsson
2 ! See http://factorcode.org/license.txt for BSD license
4 USING: combinators.short-circuit grouping kernel math
5 math.combinatorics math.functions math.order math.primes
6 math.ranges math.statistics math.vectors memoize sequences ;
14 : (sterling) ( n k -- x )
15 [ [ 1 - ] bi@ sterling ]
16 [ [ 1 - ] dip sterling ]
21 MEMO: sterling ( n k -- x )
22 2dup { [ = ] [ nip 1 = ] } 2||
23 [ 2drop 1 ] [ (sterling) ] if ;
29 : (bernoulli) ( p -- n )
31 0 [ [ nCk ] [ bernoulli * ] bi + ] with reduce
36 MEMO: bernoulli ( p -- n )
37 [ 1 ] [ (bernoulli) ] if-zero ;
39 : chi2 ( actual expected -- n )
40 0 [ dup 0 > [ [ - sq ] keep / + ] [ 2drop ] if ] 2reduce ;
45 even? [ "odd degrees of freedom" throw ] unless ;
47 : (chi2P) ( chi/2 df/2 -- p )
48 [1,b) dupd n/v cum-product swap neg e^ [ v*n sum ] keep + ;
52 : chi2P ( chi df -- p )
53 dup df-check [ 2.0 / ] [ 2 /i ] bi* (chi2P) 1.0 min ;
57 : check-jacobi ( m -- m )
58 dup { [ integer? ] [ 0 > ] [ odd? ] } 1&&
59 [ "modulus must be odd positive integer" throw ] unless ;
62 [ mod ] keep over zero? [ drop ] [
63 2dup [ sgn ] same? [ drop ] [ + ] if
69 check-jacobi [ mod' ] keep 1
73 over 8 mod' { 3 5 } member? [ neg ] when
75 2over [ 4 mod' 3 = ] both? [ neg ] when
77 ] until [ nip 1 = ] dip 0 ? ;
81 : check-legendere ( m -- m )
82 dup prime? [ "modulus must be prime positive integer" throw ] unless ;
86 : legendere ( a m -- n )
87 check-legendere jacobi ;
89 : moving-average ( seq n -- newseq )
90 <clumps> [ mean ] map ;
92 : exponential-moving-average ( seq a -- newseq )
93 [ 1 ] 2dip [ [ dupd swap - ] dip * + dup ] curry map nip ;
95 : moving-median ( u n -- v )
96 <clumps> [ median ] map ;
98 : nonzero ( seq -- seq' )
99 [ zero? not ] filter ;