Chris@49
|
1 // Copyright (C) 2009-2011 NICTA (www.nicta.com.au)
|
Chris@49
|
2 // Copyright (C) 2009-2011 Conrad Sanderson
|
Chris@49
|
3 //
|
Chris@49
|
4 // This Source Code Form is subject to the terms of the Mozilla Public
|
Chris@49
|
5 // License, v. 2.0. If a copy of the MPL was not distributed with this
|
Chris@49
|
6 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
Chris@49
|
7
|
Chris@49
|
8
|
Chris@49
|
9 //! \addtogroup running_stat
|
Chris@49
|
10 //! @{
|
Chris@49
|
11
|
Chris@49
|
12
|
Chris@49
|
13
|
Chris@49
|
14 template<typename eT>
|
Chris@49
|
15 inline
|
Chris@49
|
16 arma_counter<eT>::~arma_counter()
|
Chris@49
|
17 {
|
Chris@49
|
18 arma_extra_debug_sigprint_this(this);
|
Chris@49
|
19 }
|
Chris@49
|
20
|
Chris@49
|
21
|
Chris@49
|
22
|
Chris@49
|
23 template<typename eT>
|
Chris@49
|
24 inline
|
Chris@49
|
25 arma_counter<eT>::arma_counter()
|
Chris@49
|
26 : d_count( eT(0))
|
Chris@49
|
27 , i_count(uword(0))
|
Chris@49
|
28 {
|
Chris@49
|
29 arma_extra_debug_sigprint_this(this);
|
Chris@49
|
30 }
|
Chris@49
|
31
|
Chris@49
|
32
|
Chris@49
|
33
|
Chris@49
|
34 template<typename eT>
|
Chris@49
|
35 inline
|
Chris@49
|
36 const arma_counter<eT>&
|
Chris@49
|
37 arma_counter<eT>::operator++()
|
Chris@49
|
38 {
|
Chris@49
|
39 if(i_count < ARMA_MAX_UWORD)
|
Chris@49
|
40 {
|
Chris@49
|
41 i_count++;
|
Chris@49
|
42 }
|
Chris@49
|
43 else
|
Chris@49
|
44 {
|
Chris@49
|
45 d_count += eT(ARMA_MAX_UWORD);
|
Chris@49
|
46 i_count = 1;
|
Chris@49
|
47 }
|
Chris@49
|
48
|
Chris@49
|
49 return *this;
|
Chris@49
|
50 }
|
Chris@49
|
51
|
Chris@49
|
52
|
Chris@49
|
53
|
Chris@49
|
54 template<typename eT>
|
Chris@49
|
55 inline
|
Chris@49
|
56 void
|
Chris@49
|
57 arma_counter<eT>::operator++(int)
|
Chris@49
|
58 {
|
Chris@49
|
59 operator++();
|
Chris@49
|
60 }
|
Chris@49
|
61
|
Chris@49
|
62
|
Chris@49
|
63
|
Chris@49
|
64 template<typename eT>
|
Chris@49
|
65 inline
|
Chris@49
|
66 void
|
Chris@49
|
67 arma_counter<eT>::reset()
|
Chris@49
|
68 {
|
Chris@49
|
69 d_count = eT(0);
|
Chris@49
|
70 i_count = uword(0);
|
Chris@49
|
71 }
|
Chris@49
|
72
|
Chris@49
|
73
|
Chris@49
|
74
|
Chris@49
|
75 template<typename eT>
|
Chris@49
|
76 inline
|
Chris@49
|
77 eT
|
Chris@49
|
78 arma_counter<eT>::value() const
|
Chris@49
|
79 {
|
Chris@49
|
80 return d_count + eT(i_count);
|
Chris@49
|
81 }
|
Chris@49
|
82
|
Chris@49
|
83
|
Chris@49
|
84
|
Chris@49
|
85 template<typename eT>
|
Chris@49
|
86 inline
|
Chris@49
|
87 eT
|
Chris@49
|
88 arma_counter<eT>::value_plus_1() const
|
Chris@49
|
89 {
|
Chris@49
|
90 if(i_count < ARMA_MAX_UWORD)
|
Chris@49
|
91 {
|
Chris@49
|
92 return d_count + eT(i_count + 1);
|
Chris@49
|
93 }
|
Chris@49
|
94 else
|
Chris@49
|
95 {
|
Chris@49
|
96 return d_count + eT(ARMA_MAX_UWORD) + eT(1);
|
Chris@49
|
97 }
|
Chris@49
|
98 }
|
Chris@49
|
99
|
Chris@49
|
100
|
Chris@49
|
101
|
Chris@49
|
102 template<typename eT>
|
Chris@49
|
103 inline
|
Chris@49
|
104 eT
|
Chris@49
|
105 arma_counter<eT>::value_minus_1() const
|
Chris@49
|
106 {
|
Chris@49
|
107 if(i_count > 0)
|
Chris@49
|
108 {
|
Chris@49
|
109 return d_count + eT(i_count - 1);
|
Chris@49
|
110 }
|
Chris@49
|
111 else
|
Chris@49
|
112 {
|
Chris@49
|
113 return d_count - eT(1);
|
Chris@49
|
114 }
|
Chris@49
|
115 }
|
Chris@49
|
116
|
Chris@49
|
117
|
Chris@49
|
118
|
Chris@49
|
119 //
|
Chris@49
|
120
|
Chris@49
|
121
|
Chris@49
|
122
|
Chris@49
|
123 template<typename eT>
|
Chris@49
|
124 running_stat<eT>::~running_stat()
|
Chris@49
|
125 {
|
Chris@49
|
126 arma_extra_debug_sigprint_this(this);
|
Chris@49
|
127 }
|
Chris@49
|
128
|
Chris@49
|
129
|
Chris@49
|
130
|
Chris@49
|
131 template<typename eT>
|
Chris@49
|
132 running_stat<eT>::running_stat()
|
Chris@49
|
133 : r_mean ( eT(0))
|
Chris@49
|
134 , r_var (typename running_stat<eT>::T(0))
|
Chris@49
|
135 , min_val ( eT(0))
|
Chris@49
|
136 , max_val ( eT(0))
|
Chris@49
|
137 , min_val_norm(typename running_stat<eT>::T(0))
|
Chris@49
|
138 , max_val_norm(typename running_stat<eT>::T(0))
|
Chris@49
|
139 {
|
Chris@49
|
140 arma_extra_debug_sigprint_this(this);
|
Chris@49
|
141 }
|
Chris@49
|
142
|
Chris@49
|
143
|
Chris@49
|
144
|
Chris@49
|
145 //! update statistics to reflect new sample
|
Chris@49
|
146 template<typename eT>
|
Chris@49
|
147 inline
|
Chris@49
|
148 void
|
Chris@49
|
149 running_stat<eT>::operator() (const typename running_stat<eT>::T sample)
|
Chris@49
|
150 {
|
Chris@49
|
151 arma_extra_debug_sigprint();
|
Chris@49
|
152
|
Chris@49
|
153 if( arma_isfinite(sample) == false )
|
Chris@49
|
154 {
|
Chris@49
|
155 arma_warn(true, "running_stat: sample ignored as it is non-finite" );
|
Chris@49
|
156 return;
|
Chris@49
|
157 }
|
Chris@49
|
158
|
Chris@49
|
159 running_stat_aux::update_stats(*this, sample);
|
Chris@49
|
160 }
|
Chris@49
|
161
|
Chris@49
|
162
|
Chris@49
|
163
|
Chris@49
|
164 //! update statistics to reflect new sample (version for complex numbers)
|
Chris@49
|
165 template<typename eT>
|
Chris@49
|
166 inline
|
Chris@49
|
167 void
|
Chris@49
|
168 running_stat<eT>::operator() (const std::complex< typename running_stat<eT>::T >& sample)
|
Chris@49
|
169 {
|
Chris@49
|
170 arma_extra_debug_sigprint();
|
Chris@49
|
171
|
Chris@49
|
172 arma_type_check(( is_same_type<eT, std::complex< typename running_stat<eT>::T > >::value == false ));
|
Chris@49
|
173
|
Chris@49
|
174 if( arma_isfinite(sample) == false )
|
Chris@49
|
175 {
|
Chris@49
|
176 arma_warn(true, "running_stat: sample ignored as it is non-finite" );
|
Chris@49
|
177 return;
|
Chris@49
|
178 }
|
Chris@49
|
179
|
Chris@49
|
180 running_stat_aux::update_stats(*this, sample);
|
Chris@49
|
181 }
|
Chris@49
|
182
|
Chris@49
|
183
|
Chris@49
|
184
|
Chris@49
|
185 //! set all statistics to zero
|
Chris@49
|
186 template<typename eT>
|
Chris@49
|
187 inline
|
Chris@49
|
188 void
|
Chris@49
|
189 running_stat<eT>::reset()
|
Chris@49
|
190 {
|
Chris@49
|
191 arma_extra_debug_sigprint();
|
Chris@49
|
192
|
Chris@49
|
193 // typedef typename running_stat<eT>::T T;
|
Chris@49
|
194
|
Chris@49
|
195 counter.reset();
|
Chris@49
|
196
|
Chris@49
|
197 r_mean = eT(0);
|
Chris@49
|
198 r_var = T(0);
|
Chris@49
|
199
|
Chris@49
|
200 min_val = eT(0);
|
Chris@49
|
201 max_val = eT(0);
|
Chris@49
|
202
|
Chris@49
|
203 min_val_norm = T(0);
|
Chris@49
|
204 max_val_norm = T(0);
|
Chris@49
|
205 }
|
Chris@49
|
206
|
Chris@49
|
207
|
Chris@49
|
208
|
Chris@49
|
209 //! mean or average value
|
Chris@49
|
210 template<typename eT>
|
Chris@49
|
211 inline
|
Chris@49
|
212 eT
|
Chris@49
|
213 running_stat<eT>::mean() const
|
Chris@49
|
214 {
|
Chris@49
|
215 arma_extra_debug_sigprint();
|
Chris@49
|
216
|
Chris@49
|
217 return r_mean;
|
Chris@49
|
218 }
|
Chris@49
|
219
|
Chris@49
|
220
|
Chris@49
|
221
|
Chris@49
|
222 //! variance
|
Chris@49
|
223 template<typename eT>
|
Chris@49
|
224 inline
|
Chris@49
|
225 typename running_stat<eT>::T
|
Chris@49
|
226 running_stat<eT>::var(const uword norm_type) const
|
Chris@49
|
227 {
|
Chris@49
|
228 arma_extra_debug_sigprint();
|
Chris@49
|
229
|
Chris@49
|
230 const T N = counter.value();
|
Chris@49
|
231
|
Chris@49
|
232 if(N > T(1))
|
Chris@49
|
233 {
|
Chris@49
|
234 if(norm_type == 0)
|
Chris@49
|
235 {
|
Chris@49
|
236 return r_var;
|
Chris@49
|
237 }
|
Chris@49
|
238 else
|
Chris@49
|
239 {
|
Chris@49
|
240 const T N_minus_1 = counter.value_minus_1();
|
Chris@49
|
241 return (N_minus_1/N) * r_var;
|
Chris@49
|
242 }
|
Chris@49
|
243 }
|
Chris@49
|
244 else
|
Chris@49
|
245 {
|
Chris@49
|
246 return T(0);
|
Chris@49
|
247 }
|
Chris@49
|
248 }
|
Chris@49
|
249
|
Chris@49
|
250
|
Chris@49
|
251
|
Chris@49
|
252 //! standard deviation
|
Chris@49
|
253 template<typename eT>
|
Chris@49
|
254 inline
|
Chris@49
|
255 typename running_stat<eT>::T
|
Chris@49
|
256 running_stat<eT>::stddev(const uword norm_type) const
|
Chris@49
|
257 {
|
Chris@49
|
258 arma_extra_debug_sigprint();
|
Chris@49
|
259
|
Chris@49
|
260 return std::sqrt( (*this).var(norm_type) );
|
Chris@49
|
261 }
|
Chris@49
|
262
|
Chris@49
|
263
|
Chris@49
|
264
|
Chris@49
|
265 //! minimum value
|
Chris@49
|
266 template<typename eT>
|
Chris@49
|
267 inline
|
Chris@49
|
268 eT
|
Chris@49
|
269 running_stat<eT>::min() const
|
Chris@49
|
270 {
|
Chris@49
|
271 arma_extra_debug_sigprint();
|
Chris@49
|
272
|
Chris@49
|
273 return min_val;
|
Chris@49
|
274 }
|
Chris@49
|
275
|
Chris@49
|
276
|
Chris@49
|
277
|
Chris@49
|
278 //! maximum value
|
Chris@49
|
279 template<typename eT>
|
Chris@49
|
280 inline
|
Chris@49
|
281 eT
|
Chris@49
|
282 running_stat<eT>::max() const
|
Chris@49
|
283 {
|
Chris@49
|
284 arma_extra_debug_sigprint();
|
Chris@49
|
285
|
Chris@49
|
286 return max_val;
|
Chris@49
|
287 }
|
Chris@49
|
288
|
Chris@49
|
289
|
Chris@49
|
290
|
Chris@49
|
291 //! number of samples so far
|
Chris@49
|
292 template<typename eT>
|
Chris@49
|
293 inline
|
Chris@49
|
294 typename get_pod_type<eT>::result
|
Chris@49
|
295 running_stat<eT>::count() const
|
Chris@49
|
296 {
|
Chris@49
|
297 arma_extra_debug_sigprint();
|
Chris@49
|
298
|
Chris@49
|
299 return counter.value();
|
Chris@49
|
300 }
|
Chris@49
|
301
|
Chris@49
|
302
|
Chris@49
|
303
|
Chris@49
|
304 //! update statistics to reflect new sample
|
Chris@49
|
305 template<typename eT>
|
Chris@49
|
306 inline
|
Chris@49
|
307 void
|
Chris@49
|
308 running_stat_aux::update_stats(running_stat<eT>& x, const eT sample)
|
Chris@49
|
309 {
|
Chris@49
|
310 arma_extra_debug_sigprint();
|
Chris@49
|
311
|
Chris@49
|
312 typedef typename running_stat<eT>::T T;
|
Chris@49
|
313
|
Chris@49
|
314 const T N = x.counter.value();
|
Chris@49
|
315
|
Chris@49
|
316 if(N > T(0))
|
Chris@49
|
317 {
|
Chris@49
|
318 if(sample < x.min_val)
|
Chris@49
|
319 {
|
Chris@49
|
320 x.min_val = sample;
|
Chris@49
|
321 }
|
Chris@49
|
322
|
Chris@49
|
323 if(sample > x.max_val)
|
Chris@49
|
324 {
|
Chris@49
|
325 x.max_val = sample;
|
Chris@49
|
326 }
|
Chris@49
|
327
|
Chris@49
|
328 const T N_plus_1 = x.counter.value_plus_1();
|
Chris@49
|
329 const T N_minus_1 = x.counter.value_minus_1();
|
Chris@49
|
330
|
Chris@49
|
331 // note: variance has to be updated before the mean
|
Chris@49
|
332
|
Chris@49
|
333 const eT tmp = sample - x.r_mean;
|
Chris@49
|
334
|
Chris@49
|
335 x.r_var = N_minus_1/N * x.r_var + (tmp*tmp)/N_plus_1;
|
Chris@49
|
336
|
Chris@49
|
337 x.r_mean = x.r_mean + (sample - x.r_mean)/N_plus_1;
|
Chris@49
|
338 //x.r_mean = (N/N_plus_1)*x.r_mean + sample/N_plus_1;
|
Chris@49
|
339 //x.r_mean = (x.r_mean + sample/N) * N/N_plus_1;
|
Chris@49
|
340 }
|
Chris@49
|
341 else
|
Chris@49
|
342 {
|
Chris@49
|
343 x.r_mean = sample;
|
Chris@49
|
344 x.min_val = sample;
|
Chris@49
|
345 x.max_val = sample;
|
Chris@49
|
346
|
Chris@49
|
347 // r_var is initialised to zero
|
Chris@49
|
348 // in the constructor and reset()
|
Chris@49
|
349 }
|
Chris@49
|
350
|
Chris@49
|
351 x.counter++;
|
Chris@49
|
352 }
|
Chris@49
|
353
|
Chris@49
|
354
|
Chris@49
|
355
|
Chris@49
|
356 //! update statistics to reflect new sample (version for complex numbers)
|
Chris@49
|
357 template<typename T>
|
Chris@49
|
358 inline
|
Chris@49
|
359 void
|
Chris@49
|
360 running_stat_aux::update_stats(running_stat< std::complex<T> >& x, const T sample)
|
Chris@49
|
361 {
|
Chris@49
|
362 arma_extra_debug_sigprint();
|
Chris@49
|
363
|
Chris@49
|
364 running_stat_aux::update_stats(x, std::complex<T>(sample));
|
Chris@49
|
365 }
|
Chris@49
|
366
|
Chris@49
|
367
|
Chris@49
|
368
|
Chris@49
|
369 //! alter statistics to reflect new sample (version for complex numbers)
|
Chris@49
|
370 template<typename T>
|
Chris@49
|
371 inline
|
Chris@49
|
372 void
|
Chris@49
|
373 running_stat_aux::update_stats(running_stat< std::complex<T> >& x, const std::complex<T>& sample)
|
Chris@49
|
374 {
|
Chris@49
|
375 arma_extra_debug_sigprint();
|
Chris@49
|
376
|
Chris@49
|
377 const T sample_norm = std::norm(sample);
|
Chris@49
|
378 const T N = x.counter.value();
|
Chris@49
|
379
|
Chris@49
|
380 if(N > T(0))
|
Chris@49
|
381 {
|
Chris@49
|
382 if(sample_norm < x.min_val_norm)
|
Chris@49
|
383 {
|
Chris@49
|
384 x.min_val_norm = sample_norm;
|
Chris@49
|
385 x.min_val = sample;
|
Chris@49
|
386 }
|
Chris@49
|
387
|
Chris@49
|
388 if(sample_norm > x.max_val_norm)
|
Chris@49
|
389 {
|
Chris@49
|
390 x.max_val_norm = sample_norm;
|
Chris@49
|
391 x.max_val = sample;
|
Chris@49
|
392 }
|
Chris@49
|
393
|
Chris@49
|
394 const T N_plus_1 = x.counter.value_plus_1();
|
Chris@49
|
395 const T N_minus_1 = x.counter.value_minus_1();
|
Chris@49
|
396
|
Chris@49
|
397 x.r_var = N_minus_1/N * x.r_var + std::norm(sample - x.r_mean)/N_plus_1;
|
Chris@49
|
398
|
Chris@49
|
399 x.r_mean = x.r_mean + (sample - x.r_mean)/N_plus_1;
|
Chris@49
|
400 //x.r_mean = (N/N_plus_1)*x.r_mean + sample/N_plus_1;
|
Chris@49
|
401 //x.r_mean = (x.r_mean + sample/N) * N/N_plus_1;
|
Chris@49
|
402 }
|
Chris@49
|
403 else
|
Chris@49
|
404 {
|
Chris@49
|
405 x.r_mean = sample;
|
Chris@49
|
406 x.min_val = sample;
|
Chris@49
|
407 x.max_val = sample;
|
Chris@49
|
408 x.min_val_norm = sample_norm;
|
Chris@49
|
409 x.max_val_norm = sample_norm;
|
Chris@49
|
410
|
Chris@49
|
411 // r_var is initialised to zero
|
Chris@49
|
412 // in the constructor and reset()
|
Chris@49
|
413 }
|
Chris@49
|
414
|
Chris@49
|
415 x.counter++;
|
Chris@49
|
416 }
|
Chris@49
|
417
|
Chris@49
|
418
|
Chris@49
|
419
|
Chris@49
|
420 //! @}
|