/
stats.4th
179 lines (142 loc) · 4.52 KB
/
stats.4th
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
\ stats.4th
\
\ Compute the mean and variance of a set of floating point numbers
\
\ Examples: See test code
\
\ Requires:
\ fsl-util.4th
\
\ Copyright (c) 1998, 2003 Krishna Myneni
\ Revisions:
\ 1998-12-22
\ 2003-11-06 km; added probability density function calculations:
\ PGAUSS, AGAUSS, AREA_GAUSS
\ 2007-11-09 km; revised for FSL compatibility
\ 2010-10-26 km; renamed }mean and }variance to }fmean and
\ }fvariance for consistency with array-utils
\ naming convention; also use Private: and
\ Public: to hide transient variables
\ 2011-09-16 km; use Neal Bridges' anonymous modules.
\ 2012-02-19 km; use KM/DNW's modules library
\ 2013-03-12 km; fix stack comment for .STATS
\ 2017-04-21 km; added }IMEAN and }IVARIANCE ; made .stats
\ a private word and provided }ISTATS. and }FSTATS.
\ 2023-12-05 km; replaced FDUP F* with FSQUARE
CR .( STATS V1.1h 05 December 2023 KM )
BEGIN-MODULE
BASE @
DECIMAL
Private:
0 ptr arr{
fvariable mu
fvariable sigma2
: .stats ( -- )
mu f@ ." Mean = " F. cr
sigma2 f@ ." Variance = " F. cr
;
Public:
0 value Nsamples
: }imean ( u 'a -- fmu | mean of integer array )
TO arr{ TO Nsamples
0 Nsamples 0 ?DO arr{ I } @ + LOOP s>f Nsamples s>f F/
fdup mu f! ;
: }ivariance ( u 'a -- fsigma2 | variance of integer array )
}imean
0e Nsamples 0 ?DO fover arr{ I } @ s>f F- fsquare F+ LOOP
fswap fdrop
Nsamples 1- s>f F/
fdup sigma2 F! ;
: }fmean ( u 'a -- fmu | mean of floating point array )
TO arr{ TO Nsamples
0e Nsamples 0 ?DO arr{ I } F@ F+ LOOP Nsamples s>f F/
fdup mu f! ;
: }fvariance ( u 'a -- fsigma2 | variance of floating point array )
}fmean
0e Nsamples 0 ?DO fover arr{ I } F@ F- fsquare F+ LOOP
fswap fdrop
Nsamples 1- s>f F/
fdup sigma2 F! ;
: }istats. ( u 'a -- | compute and print statistics for integer array )
}ivariance fdrop .stats ;
: }fstats. ( u 'a -- | compute and print statistics for fp array )
}fvariance fdrop .stats ;
2e fsqrt fconstant sqrt{2}
\ The words pgauss and agauss are translated from corresponding
\ Fortran functions in P.R. Bevington, Data Reduction and Error
\ Analysis for the Physical Sciences, 1969, McGraw-Hill.
: pgauss ( fx fmu fsigma -- fpdf | evaluate gaussian probability density at fx)
fdup F0= ABORT" PGAUSS: zero sigma value"
2>r F- 2r@ F/ fsquare 2e F/ fnegate fexp
0.3989422804e F* 2r> F/ ;
\ Evaluate area between the limits (mu - z*sigma) to (mu + z*sigma)
\ where z = |x-mu|/sigma, e.g. if mu = 0, the area is computed
\ between -x and +x.
Private:
fvariable agauss_term
fvariable agauss_sum
fvariable agauss_denom
Public:
: agauss ( fx fmu fsigma -- farea )
fdup F0= ABORT" AGAUSS: zero sigma value"
2>r F- 2r> F/ fabs \ -- z
fdup 0e F<= IF fdrop 0e exit THEN \ -- z
fdup sqrt{2} F/ \ -- z term
fdup agauss_term F! agauss_sum F! \ -- z
fsquare 2e F/ \ -- y2
\ Accumulate sum of terms
1e agauss_denom F!
BEGIN
fdup 2e F* \ -- y2 2*y2
agauss_denom F@ 2e F+ fdup agauss_denom F! F/ \ -- y2 2*y2/denom
agauss_term F@ F* fdup
agauss_sum F@ F+ agauss_sum F!
fdup agauss_term F!
agauss_sum F@ F/ 1e-10 F<=
UNTIL
fnegate fexp agauss_sum F@ F* 1.128379167e F*
;
\ Return area between x1 and x2 for Gaussian probability density function
Private:
fvariable agauss_mu
fvariable agauss_sigma
Public:
: area_gauss ( fx1 fx2 fmu fsigma -- farea )
fdup F0= ABORT" AREA_GAUSS: zero sigma value"
agauss_sigma F!
agauss_mu F! \ -- fx1 fx2
agauss_mu F@ F- fswap
agauss_mu F@ F- fswap \ -- fx1-fmu fx2-fmu
fover fover F* 0e F> >r
fabs 0e agauss_sigma F@ agauss fswap
fabs 0e agauss_sigma F@ agauss
r>
IF F- fabs \ (x1-mu) and (x2-mu) have same sign?
ELSE F+ \ (x1-mu) and (x2-mu) have opposite sign
THEN 2e F/
;
BASE !
END-MODULE
TEST-CODE? [IF] \ test code ==============================================
[undefined] T{ [IF] include ttester.4th [THEN]
BASE @ DECIMAL
1e-7 rel-near F!
1e-7 abs-near F!
set-near
(
\ Test code for Gaussian probability functions
: test_gauss \ -- | print a table of values for mu = 0, sigma =1
0e
BEGIN
cr
fdup f. 2 spaces
fdup 0e 1e pgauss f. 2 spaces
fdup 0e 1e agauss f. 2 spaces
fdup fnegate fover 0e 1e area_gauss f.
0.1e f+
fdup 3e f>
UNTIL
fdrop cr ;
)
BASE !
[THEN]