Home My Page Projects Code Snippets Project Openings SML/NJ
Summary Activity Forums Tracker Lists Tasks Docs Surveys News SCM Files

SCM Repository

[smlnj] Annotation of /sml/trunk/src/smlnj-lib/Util/univariate-stats.sml
ViewVC logotype

Annotation of /sml/trunk/src/smlnj-lib/Util/univariate-stats.sml

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1720 - (view) (download)

1 : mblume 1720 (* univariate-stats.sml
2 :     *
3 :     * Some statistical functions on unweighted univariate samples.
4 :     *
5 :     * Copyright (c) 2004 by The Fellowship of SML/NJ
6 :     *
7 :     * Author: Matthias Blume (blume@tti-c.org)
8 :     *)
9 :     structure UnivariateStats :> sig
10 :    
11 :     (* We distinguish between two kinds of samples. Only the "heavy"
12 :     * kind permits calculation of average deviation and median.
13 :     * It is also considerably more expensive because it keeps an
14 :     * array of all points while the "light" variety is constant-size. *)
15 :     type light type heavy
16 :    
17 :     type 'a sample (* light or heavy *)
18 :     type 'a evaluation (* light or heavy *)
19 :    
20 :     (* Samples are built incrementally by adding points to an initially
21 :     * empty sample: *)
22 :     val lempty : light sample
23 :     val hempty : unit -> heavy sample
24 :     val ladd : real * light sample -> light sample (* constant *)
25 :     val hadd : real * heavy sample -> heavy sample (* amortized constant *)
26 :    
27 :     (* Evaluate the sample; this completes all the expensive work except
28 :     * for things that depend on "heavy" samples: *)
29 :     val evaluate : 'a sample -> 'a evaluation (* constant *)
30 :    
31 :     (* extracting of "cheap" information (constant-time): *)
32 :     val N : 'a evaluation -> int
33 :     val n : 'a evaluation -> real (* N as real *)
34 :     val mean : 'a evaluation -> real
35 :     val variance : 'a evaluation -> real
36 :     val standardDeviation : 'a evaluation -> real
37 :     val skew : 'a evaluation -> real
38 :     val kurtosis : 'a evaluation -> real
39 :    
40 :     (* extracting of "expensive" information: *)
41 :     val median : heavy evaluation -> real (* randomized linear *)
42 :     val averageDeviation : heavy evaluation -> real (* linear *)
43 :    
44 :     end = struct
45 :    
46 :     infix 8 $ val op $ = Unsafe.Array.sub
47 :     infix 3 <- fun (a, i) <- x = Unsafe.Array.update (a, i, x)
48 :    
49 :     type light = unit
50 :     type heavy = real array * int
51 :     type 'a sample = ('a * int * real * real * real * real)
52 :     type 'a evaluation = ('a * int * real * real * real * real * real * real)
53 :    
54 :     fun insert (x, n, (a, sz)) =
55 :     let val (a, sz) =
56 :     if n<sz then (a, sz)
57 :     else let val sz = sz+sz
58 :     val b=Array.tabulate(sz,fn i=>if i<n then a$i else 0.0)
59 :     in (b, sz) end
60 :     in (a,n)<-x; (a, sz) end
61 :    
62 :     val SZ = 1024 (* minimum allocated size of heavy array *)
63 :     val lempty = ((), 0, 0.0, 0.0, 0.0, 0.0)
64 :     fun hempty () = ((Array.array (SZ, 0.0), SZ), 0, 0.0, 0.0, 0.0, 0.0)
65 :    
66 :     fun ladd (x:real, ((), n, sx4, sx3, sx2, sx1)) =
67 :     let val x2 = x*x val (x3, x4) = (x2*x, x2*x2)
68 :     in ((), n+1, sx4+x4, sx3+x3, sx2+x2, sx1+x) end
69 :    
70 :     fun hadd (x:real, ((a, sz), n, sx4, sx3, sx2, sx1)) =
71 :     let val x2 = x*x val (x3, x4) = (x2*x, x2*x2)
72 :     val (a, sz) =
73 :     if n < sz then (a, sz)
74 :     else let val sz = sz+sz
75 :     val b = Array.tabulate
76 :     (sz, fn i => if i<n then a$i else 0.0)
77 :     in (b, sz) end
78 :     in (a,n)<-x;
79 :     ((a, sz), n+1, sx4+x4, sx3+x3, sx2+x2, sx1+x)
80 :     end
81 :    
82 :     fun evaluate (state, ni, sx4, sx3, sx2, sx1) =
83 :     let val n = real ni val n' = n - 1.0
84 :     val m = sx1/n val m2 = m*m val m3 = m2*m
85 :     val sd2 = (sx2 + m*(n*m-2.0*sx1))/n'
86 :     val sd = Math.sqrt sd2 val (sd3, sd4) = (sd*sd2, sd2*sd2)
87 :     val sk = (sx3-m*(3.0*(sx2-sx1*m)+n*m2))/(n*sd3)
88 :     val k = ((sx4+m*(6.0*sx2*m-4.0*(sx3+sx1*m2)+n*m3))/(n*sd4))-3.0
89 :     in (state, ni, n, m, sd2, sd, sk, k) end
90 :    
91 :     fun N (state, ni, nr, m, sd2, sd, sk, k) = ni
92 :     fun n (state, ni, nr, m, sd2, sd, sk, k) = nr
93 :     fun mean (state, ni, nr, m, sd2, sd, sk, k) = m
94 :     fun variance (state, ni, nr, m, sd2, sd, sk, k) = sd2
95 :     fun standardDeviation (state, ni, nr, m, sd2, sd, sk, k) = sd
96 :     fun skew (state, ni, nr, m, sd2, sd, sk, k) = sk
97 :     fun kurtosis (state, ni, nr, m, sd2, sd, sk, k) = k
98 :    
99 :     fun median ((a, sz), ni, nr, m, sd2, sd, sk, k) =
100 :     RealOrderStats.median' (ArraySlice.slice (a, 0, SOME ni))
101 :    
102 :     fun averageDeviation ((a, sz), ni, nr, m, sd2, sd, sk, k) =
103 :     let fun ad (i, ds) = if i>ni then ds/nr else ad (i+1, ds + abs(a$i-m))
104 :     in ad (0, 0.0) end
105 :     end

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0