2
2
//
3
3
// FILE: Gauss.h
4
4
// AUTHOR: Rob Tillaart
5
- // VERSION: 0.1.1
5
+ // VERSION: 0.2.0
6
6
// PURPOSE: Library for the Gauss probability math.
7
7
// DATE: 2023-07-06
8
8
9
9
10
10
#include " Arduino.h"
11
- #include " MultiMap.h"
12
11
13
- #define GAUSS_LIB_VERSION (F(" 0.1.1 " ))
12
+ #define GAUSS_LIB_VERSION (F(" 0.2.0 " ))
14
13
15
14
16
15
class Gauss
@@ -19,16 +18,15 @@ class Gauss
19
18
Gauss ()
20
19
{
21
20
_mean = 0 ;
22
- _stddev = 1 ;
23
21
_reciprokeSD = 1 ;
24
22
}
25
23
26
24
// stddev should be positive.
27
25
bool begin (float mean = 0 , float stddev = 1 )
28
26
{
29
27
_mean = mean;
30
- _stddev = stddev; // should be positive
31
- _reciprokeSD = 1.0 / _stddev ;
28
+ if ( stddev == 0 ) _reciprokeSD = NAN;
29
+ else _reciprokeSD = 1.0 / stddev ;
32
30
return (stddev > 0 );
33
31
}
34
32
@@ -41,14 +39,13 @@ class Gauss
41
39
42
40
float getStdDev ()
43
41
{
44
- return _stddev ;
42
+ return 1.0 / _reciprokeSD ;
45
43
}
46
44
47
45
48
46
float P_smaller (float value)
49
47
{
50
- if (_stddev == 0 ) return NAN;
51
- // normalize(value)
48
+ if (_reciprokeSD == NAN) return NAN;
52
49
return _P_smaller ((value - _mean) * _reciprokeSD);
53
50
}
54
51
@@ -61,18 +58,25 @@ class Gauss
61
58
62
59
float P_between (float p, float q)
63
60
{
64
- if (_stddev == 0 ) return NAN;
61
+ if (_reciprokeSD == NAN ) return NAN;
65
62
if (p >= q) return 0 ;
66
63
return P_smaller (q) - P_smaller (p);
67
64
}
68
65
69
66
67
+ float P_outside (float p, float q)
68
+ {
69
+ return 1.0 - P_between (p, q);
70
+ }
71
+
72
+
70
73
float P_equal (float value)
71
74
{
72
- if (_stddev == 0 ) return NAN;
75
+ if (_reciprokeSD == NAN ) return NAN;
73
76
float n = (value - _mean) * _reciprokeSD;
77
+ // gain of ~10% if we allocate a global var for 'constant' c
74
78
float c = _reciprokeSD * (1.0 / sqrt (TWO_PI));
75
- return c * exp (-0.5 * n * n);
79
+ return c * exp (-0.5 * ( n * n) );
76
80
}
77
81
78
82
@@ -88,18 +92,32 @@ class Gauss
88
92
}
89
93
90
94
95
+ float denormalize (float value)
96
+ {
97
+ return value / _reciprokeSD + _mean;
98
+ }
99
+
100
+
91
101
float bellCurve (float value)
92
102
{
93
103
return P_equal (value);
94
104
}
95
105
96
106
107
+ float CDF (float value)
108
+ {
109
+ return P_smaller (value);
110
+ }
111
+
112
+
97
113
98
114
private:
99
115
100
116
float _P_smaller (float x)
101
117
{
102
118
// NORM.DIST(mean, stddev, x, true)
119
+ // these points correspond with
120
+ // 0.0 .. 3.0 in steps of 0.1 followed by 4.0, 5.0 and 6.0
103
121
float __gauss[] = {
104
122
0.50000000 , 0.53982784 , 0.57925971 , 0.61791142 ,
105
123
0.65542174 , 0.69146246 , 0.72574688 , 0.75803635 ,
@@ -112,24 +130,37 @@ class Gauss
112
130
0.99999971 , 1.00000000
113
131
};
114
132
115
- // 0..60000 uint16_t = 68 bytes less
116
- float __z[] = {
117
- 0.0 , 0.1 , 0.2 , 0.3 , 0.4 , 0.5 , 0.6 , 0.7 ,
118
- 0.8 , 0.9 , 1.0 , 1.1 , 1.2 , 1.3 , 1.4 , 1.5 ,
119
- 1.6 , 1.7 , 1.8 , 1.9 , 2.0 , 2.1 , 2.2 , 2.3 ,
120
- 2.4 , 2.5 , 2.6 , 2.7 , 2.8 , 2.9 , 3.0 , 4 ,0 ,
121
- 5.0 , 6.0
122
- };
123
-
124
- // a dedicated MultiMap could exploit the fact that
125
- // the __z[] array is largely equidistant.
126
- // that could remove the __z[] array (almost) completely.
127
- if (x < 0 ) return 1.0 - multiMap<float >(-x, __z, __gauss, 34 );
128
- return multiMap<float >(x, __z, __gauss, 34 );
133
+ bool neg = false ;
134
+ if (x < 0 )
135
+ {
136
+ neg = true ;
137
+ x = -x;
138
+ }
139
+
140
+ if (x >= 6.0 )
141
+ {
142
+ if (neg) return 0.0 ;
143
+ return 1.0 ;
144
+ }
145
+
146
+ if (x <= 3.0 )
147
+ {
148
+ int idx = x * 10 ;
149
+ float rv = __gauss[idx] + ((x * 10 ) - idx) * (__gauss[idx+1 ] - __gauss[idx]);
150
+ if (neg) return 1.0 - rv;
151
+ return rv;
152
+ }
153
+
154
+ // 3.0 .. 6.0
155
+ int xint = x;
156
+ int idx = 27 + xint;
157
+ float rv = __gauss[idx] + (x - xint) * (__gauss[idx+1 ] - __gauss[idx]);
158
+ if (neg) return 1.0 - rv;
159
+ return rv;
129
160
}
130
161
131
162
float _mean = 0 ;
132
- float _stddev = 1 ; // not needed as _reciprokeSD holds same info?
163
+ // reciprokeSD = 1.0 / stddev is faster in most math (MUL vs DIV)
133
164
float _reciprokeSD = 1 ;
134
165
};
135
166
0 commit comments