22//
33// FILE: Gauss.h
44// AUTHOR: Rob Tillaart
5- // VERSION: 0.1.1
5+ // VERSION: 0.2.0
66// PURPOSE: Library for the Gauss probability math.
77// DATE: 2023-07-06
88
99
1010#include " Arduino.h"
11- #include " MultiMap.h"
1211
13- #define GAUSS_LIB_VERSION (F(" 0.1.1 " ))
12+ #define GAUSS_LIB_VERSION (F(" 0.2.0 " ))
1413
1514
1615class Gauss
@@ -19,16 +18,15 @@ class Gauss
1918 Gauss ()
2019 {
2120 _mean = 0 ;
22- _stddev = 1 ;
2321 _reciprokeSD = 1 ;
2422 }
2523
2624 // stddev should be positive.
2725 bool begin (float mean = 0 , float stddev = 1 )
2826 {
2927 _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 ;
3230 return (stddev > 0 );
3331 }
3432
@@ -41,14 +39,13 @@ class Gauss
4139
4240 float getStdDev ()
4341 {
44- return _stddev ;
42+ return 1.0 / _reciprokeSD ;
4543 }
4644
4745
4846 float P_smaller (float value)
4947 {
50- if (_stddev == 0 ) return NAN;
51- // normalize(value)
48+ if (_reciprokeSD == NAN) return NAN;
5249 return _P_smaller ((value - _mean) * _reciprokeSD);
5350 }
5451
@@ -61,18 +58,25 @@ class Gauss
6158
6259 float P_between (float p, float q)
6360 {
64- if (_stddev == 0 ) return NAN;
61+ if (_reciprokeSD == NAN ) return NAN;
6562 if (p >= q) return 0 ;
6663 return P_smaller (q) - P_smaller (p);
6764 }
6865
6966
67+ float P_outside (float p, float q)
68+ {
69+ return 1.0 - P_between (p, q);
70+ }
71+
72+
7073 float P_equal (float value)
7174 {
72- if (_stddev == 0 ) return NAN;
75+ if (_reciprokeSD == NAN ) return NAN;
7376 float n = (value - _mean) * _reciprokeSD;
77+ // gain of ~10% if we allocate a global var for 'constant' c
7478 float c = _reciprokeSD * (1.0 / sqrt (TWO_PI));
75- return c * exp (-0.5 * n * n);
79+ return c * exp (-0.5 * ( n * n) );
7680 }
7781
7882
@@ -88,18 +92,32 @@ class Gauss
8892 }
8993
9094
95+ float denormalize (float value)
96+ {
97+ return value / _reciprokeSD + _mean;
98+ }
99+
100+
91101 float bellCurve (float value)
92102 {
93103 return P_equal (value);
94104 }
95105
96106
107+ float CDF (float value)
108+ {
109+ return P_smaller (value);
110+ }
111+
112+
97113
98114private:
99115
100116 float _P_smaller (float x)
101117 {
102118 // 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
103121 float __gauss[] = {
104122 0.50000000 , 0.53982784 , 0.57925971 , 0.61791142 ,
105123 0.65542174 , 0.69146246 , 0.72574688 , 0.75803635 ,
@@ -112,24 +130,37 @@ class Gauss
112130 0.99999971 , 1.00000000
113131 };
114132
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;
129160 }
130161
131162 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)
133164 float _reciprokeSD = 1 ;
134165};
135166
0 commit comments