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
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
|
#include "c23-fallback.h"
#include <tgmath.h>
#include <stdio.h>
/**
** @file
** @brief Collect some simple statistics online as we go.
**
** This uses a generalization of Welford's trick to compute running
** mean and variance. See
**
** Philippe Pébay. Formulas for robust, one-pass parallel computation
** of covariances and arbitrary-order statistical moments. Technical
** Report SAND2008-6212, SANDIA, 2008. URL
** http://prod.sandia.gov/techlib/access-control.cgi/2008/086212.pdf.
**/
/**
** @brief A simple data structure to collect the 0th to 3rd moment of
** a statistic.
**
** @warning Since this also uses a @c double for the number of
** samples, the validity of all this is restricted to about
** @f$2^{50} \approx 10^{15}@f$ samples.
**/
struct stats {
double moment[4];
};
typedef struct stats stats;
/**
** @brief Return the number of samples that had been entered into the
** statistic @a c.
**/
inline
double stats_samples(stats c[static 1]) [[__unsequenced__]] {
return c->moment[0];
}
/**
** @brief Return the mean value of the samples that had been entered
** into the statistic @a c.
**/
inline
double stats_mean(stats c[static 1]) [[__unsequenced__]] {
return c->moment[1];
}
/**
** @brief Return the variance of the samples that had been entered
** into the statistic @a c.
**/
inline
double stats_var(stats c[static 1]) [[__reproducible__]] {
return c->moment[2]/stats_samples(c);
}
// Clang has this notorious design flaw in <tgmath.h> where they
// declare interfaces as being static. In some cases this results in
// an avalanche of confusing messages of static functions that would
// be used in an inline function.
#if __clang_major__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wstatic-in-inline"
#endif
/**
** @brief Return the standard deviation of the samples that had been
** entered into the statistic @a c.
**/
inline
double stats_sdev(stats c[static 1]) [[__reproducible__]] {
return sqrt(stats_var(c));
}
/**
** @brief Return the relative standard deviation of the samples that
** had been entered into the statistic @a c.
**/
inline
double stats_rsdev(stats c[static 1]) [[__reproducible__]] {
return sqrt(stats_var(c))/stats_mean(c);
}
/**
** @brief Return the normalized skew of the samples that had been
** entered into the statistic @a c.
**/
inline
double stats_skew(stats c[static 1]) [[__reproducible__]] {
double var = stats_var(c);
return (c->moment[3]/pow(var, 1.5))/stats_samples(c);
}
/**
** @brief Return the unbiased variance of the samples that had been
** entered into the statistic @a c.
**
** Use Bessel's correction to have an estimation of the unbiased
** variance of the overall population.
**/
inline
double stats_var_unbiased(stats c[static 1]) [[__reproducible__]] {
return c->moment[2]/(stats_samples(c)-1);
}
/**
** @brief Return the unbiased standard deviation of the samples that
** had been entered into the statistic @a c.
**
** Use Bessel's correction to have an less biased estimation of the
** variance of the overall population.
**/
inline
double stats_sdev_unbiased(stats c[static 1]) [[__reproducible__]] {
return sqrt(stats_var_unbiased(c));
}
#if __clang_major__
#pragma clang diagnostic pop
#endif
/**
** @brief Return the unbiased relative standard deviation of the
** samples that had been entered into the statistic @a c.
**/
inline
double stats_rsdev_unbiased(stats c[static 1]) [[__reproducible__]] {
return stats_rsdev(c)*(1+1/(4*stats_samples(c)));
}
/**
** @brief Add value @a val to the statistic @a c.
**
** @c moments is the number of statistic moments that is collected,
** it has to be between `0` and `3`, including.
**/
inline
void stats_collect(stats c[static 1],
double val, unsigned moments)
[[__reproducible__]] {
double n = stats_samples(c);
double n0 = n-1;
double n1 = n+1;
double delta0 = 1;
double delta = val - stats_mean(c);
double delta1 = delta/n1;
double delta2 = delta1*delta*n;
switch (moments) {
default:
c->moment[3] += (delta2*n0 - 3*c->moment[2])*delta1;
[[__fallthrough__]];
case 2:
c->moment[2] += delta2;
[[__fallthrough__]];
case 1:
c->moment[1] += delta1;
[[__fallthrough__]];
case 0:
c->moment[0] += delta0;
}
}
/**
** @brief Add value @a val to the statistic @a c.
**
** Only the number of samples is collected.
**/
inline
void stats_collect0(stats c[static 1],
double val)
[[__reproducible__]] {
stats_collect(c, val, 0);
}
/**
** @brief Add value @a val to the statistic @a c.
**
** Only the number of samples and mean is collected.
**/
inline
void stats_collect1(stats c[static 1],
double val)
[[__reproducible__]] {
stats_collect(c, val, 1);
}
/**
** @brief Add value @a val to the statistic @a c.
**
** The number of samples, mean and standard deviation is collected.
**/
inline
void stats_collect2(stats c[static 1],
double val)
[[__reproducible__]] {
stats_collect(c, val, 2);
}
/**
** @brief Add value @a val to the statistic @a c.
**
** All statistics up to the skew are collected.
**/
inline
void stats_collect3(stats c[static 1],
double val)
[[__reproducible__]] {
stats_collect(c, val, 3);
}
|
