basis/math/statistics/statistics.factor

   1 ! Copyright (C) 2008 Doug Coleman, Michael Judge.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: arrays combinators kernel math math.analysis
   4 math.functions math.order sequences sorting ;
   5 IN: math.statistics
   6
   7 : mean ( seq -- n )
   8     [ sum ] [ length ] bi / ;
   9
  10 : geometric-mean ( seq -- n )
  11     [ length ] [ product ] bi nth-root ;
  12
  13 : harmonic-mean ( seq -- n )
  14     [ recip ] sigma recip ;
  15
  16 : median ( seq -- n )
  17     natural-sort dup length even? [
  18         [ midpoint@ dup 1- 2array ] keep nths mean
  19     ] [
  20         [ midpoint@ ] keep nth
  21     ] if ;
  22
  23 : minmax ( seq -- min max )
  24     #! find the min and max of a seq in one pass
  25     [ 1/0. -1/0. ] dip [ tuck [ min ] [ max ] 2bi* ] each ;
  26
  27 : range ( seq -- n )
  28     minmax swap - ;
  29
  30 : var ( seq -- x )
  31     #! normalize by N-1
  32     dup length 1 <= [
  33         drop 0
  34     ] [
  35         [ [ mean ] keep [ - sq ] with sigma ] keep
  36         length 1- /
  37     ] if ;
  38
  39 : std ( seq -- x )
  40     var sqrt ;
  41
  42 : ste ( seq -- x )
  43     [ std ] [ length ] bi sqrt / ;
  44
  45 : ((r)) ( mean(x) mean(y) {x} {y} -- (r) )
  46     ! finds sigma((xi-mean(x))(yi-mean(y))
  47     0 [ [ [ pick ] dip swap - ] bi@ * + ] 2reduce 2nip ;
  48
  49 : (r) ( mean(x) mean(y) {x} {y} sx sy -- r )
  50     * recip [ [ ((r)) ] keep length 1- / ] dip * ;
  51
  52 : [r] ( {{x,y}...} -- mean(x) mean(y) {x} {y} sx sy )
  53     first2 [ [ [ mean ] bi@ ] 2keep ] 2keep [ std ] bi@ ;
  54
  55 : r ( {{x,y}...} -- r )
  56     [r] (r) ;
  57
  58 : r^2 ( {{x,y}...} -- r )
  59     r sq ;
  60
  61 : least-squares ( {{x,y}...} -- alpha beta )
  62     [r] { [ 2dup ] [ ] [ ] [ ] [ ] } spread
  63     ! stack is mean(x) mean(y) mean(x) mean(y) {x} {y} sx sy
  64     [ (r) ] 2keep ! stack is mean(x) mean(y) r sx sy
  65     swap / * ! stack is mean(x) mean(y) beta
  66     [ swapd * - ] keep ;
  67