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 math.functions sequences
   4     sequences.lib 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 : range ( seq -- n )
  24     minmax swap - ;
  25
  26 : var ( seq -- x )
  27     #! normalize by N-1
  28     dup length 1 <= [
  29         drop 0
  30     ] [
  31         [ [ mean ] keep [ - sq ] with sigma ] keep
  32         length 1- /
  33     ] if ;
  34
  35 : std ( seq -- x )
  36     var sqrt ;
  37
  38 : ste ( seq -- x )
  39     [ std ] [ length ] bi sqrt / ;
  40
  41 : ((r)) ( mean(x) mean(y) {x} {y} -- (r) )
  42     ! finds sigma((xi-mean(x))(yi-mean(y))
  43     0 [ [ [ pick ] dip swap - ] bi@ * + ] 2reduce 2nip ;
  44
  45 : (r) ( mean(x) mean(y) {x} {y} sx sy -- r )
  46     * recip [ [ ((r)) ] keep length 1- / ] dip * ;
  47
  48 : [r] ( {{x,y}...} -- mean(x) mean(y) {x} {y} sx sy )
  49     first2 [ [ [ mean ] bi@ ] 2keep ] 2keep [ std ] bi@ ;
  50
  51 : r ( {{x,y}...} -- r )
  52     [r] (r) ;
  53
  54 : r^2 ( {{x,y}...} -- r )
  55     r sq ;
  56
  57 : least-squares ( {{x,y}...} -- alpha beta )
  58     [r] { [ 2dup ] [ ] [ ] [ ] [ ] } spread
  59     ! stack is mean(x) mean(y) mean(x) mean(y) {x} {y} sx sy
  60     [ (r) ] 2keep ! stack is mean(x) mean(y) r sx sy
  61     swap / * ! stack is mean(x) mean(y) beta
  62     [ swapd * - ] keep ;
  63