vec_stats.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3
4  Copyright (C) 2006-2017, Andrew W. Steiner
5
6  This file is part of O2scl.
7
8  O2scl is free software; you can redistribute it and/or modify
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20
21  -------------------------------------------------------------------
22 */
23 #ifndef O2SCL_VEC_STATS_H
24 #define O2SCL_VEC_STATS_H
25
26 /** \file vec_stats.h
27  \brief Statistical functions for vector types
28
29  This file contains several function templates for computing
30  statistics of vectors of double-precision data. It includes mean,
31  median, variance, standard deviation, covariance, correlation, and
32  other functions.
33
34  No additional range checking is done on the vectors.
35
36  \future Consider generalizing to other data types.
37 */
38
39 #include <o2scl/err_hnd.h>
40 #include <o2scl/vector.h>
41
42 #ifndef DOXYGEN_NO_O2NS
43 namespace o2scl {
44 #endif
45
46  /// \name Vector functions
47  //@{
48  /** \brief Compute the mean of the first \c n elements of a vector
49
50  This function produces the same results
51  as <tt>gsl_stats_mean()</tt>.
52
53  If \c n is zero, this function will return zero.
54  */
55  template<class vec_t> double vector_mean(size_t n, const vec_t &data) {
56  long double mean=0.0;
57  for(size_t i=0;i<n;i++) {
58  mean+=(data[i]-mean)/(i+1);
59  }
60  return mean;
61  }
62
63  /** \brief Compute the mean of all of the vector elements
64
65  This function uses <tt>size()</tt> to determine the vector size
66  and produces the same results as <tt>gsl_stats_mean()</tt>.
67
68  If the vector size is zero, this function will return zero.
69  */
70  template<class vec_t> double vector_mean(const vec_t &data) {
71  return vector_mean(data.size(),data);
72  }
73
74  /** \brief Compute variance with specified mean known in advance
75
76  This function computes
77  \f[
78  \frac{1}{N} \sum_{i} \left( x_i - \mu \right)^2
79  \f]
80  where the value of \f$\mu \f$ is given in \c mean.
81
82  This function produces the same results as
83  <tt>gsl_stats_variance_with_fixed_mean()</tt>. IF \c n is zero,
84  this function returns zero.
85  */
86  template<class vec_t>
87  double vector_variance_fmean(size_t n, const vec_t &data, double mean) {
88  long double var=0.0;
89  for(size_t i=0;i<n;i++) {
90  long double delta=(data[i]-mean);
91  var+=(delta*delta-var)/(i+1);
92  }
93  return var;
94  }
95
96  /** \brief Compute variance with specified mean known in advance
97
98  This function computes
99  \f[
100  \frac{1}{N} \sum_{i} \left( x_i - \mu \right)^2
101  \f]
102  where the value of \f$\mu \f$ is given in \c mean.
103
104  This function produces the same results as
105  <tt>gsl_stats_variance_with_fixed_mean()</tt>. If the vector
106  size is zero, this function will return zero.
107  */
108  template<class vec_t>
109  double vector_variance_fmean(const vec_t &data, double mean) {
110  return vector_variance_fmean(data.size(),data,mean);
111  }
112
113  /** \brief Compute the variance with specified mean
114
115  This function computes
116  \f[
117  \frac{1}{N-1} \sum_{i} \left( x_i - \mu \right)^2
118  \f]
119  where the value of \f$\mu \f$ is given in \c mean.
120
121  This function produces the same results
122  as <tt>gsl_stats_variance_m</tt>.
123
124  If \c n is 0 or 1, this function will call the error
125  handler.
126  */
127  template<class vec_t>
128  double vector_variance(size_t n, const vec_t &data, double mean) {
129
130  if (n<2) {
131  O2SCL_ERR2("Cannot compute variance with less than 2 elements",
132  " in vector_variance().",exc_einval);
133  }
134
135  double var=vector_variance_fmean<vec_t>(n,data,mean);
136  return var*n/(n-1);
137  }
138
139  /** \brief Compute the variance with specified mean
140
141  This function computes
142  \f[
143  \frac{1}{N-1} \sum_{i} \left( x_i - \mu \right)^2
144  \f]
145  where the value of \f$\mu \f$ is given in \c mean.
146
147  This function produces the same results
148  as <tt>gsl_stats_variance_m</tt>.
149
150  If \c n is 0 or 1, this function will call the error
151  handler.
152  */
153  template<class vec_t>
154  double vector_variance(const vec_t &data, double mean) {
155  return vector_variance(data.size(),data,mean);
156  }
157
158  /** \brief Compute the variance
159
160  This function computes
161  \f[
162  \frac{1}{N-1} \sum_{i} \left( x_i - \mu \right)^2
163  \f]
164  where \f$\mu \f$ is the mean computed with \ref vector_mean().
165
166  This function produces the same results
167  as <tt>gsl_stats_variance</tt>.
168
169  If \c n is 0 or 1, this function will call the error handler.
170  */
171  template<class vec_t> double vector_variance(size_t n, const vec_t &data) {
172
173  if (n<2) {
174  O2SCL_ERR2("Cannot compute variance with less than 2 elements",
175  " in vector_variance().",exc_einval);
176  }
177
178  double mean=vector_mean<vec_t>(n,data);
179  double var=vector_variance_fmean<vec_t>(n,data,mean);
180  return var*n/(n-1);
181  }
182
183  /** \brief Compute the variance
184
185  This function computes
186  \f[
187  \frac{1}{N-1} \sum_{i} \left( x_i - \mu \right)^2
188  \f]
189  where \f$\mu \f$ is the mean computed with \ref vector_mean().
190
191  This function produces the same results
192  as <tt>gsl_stats_variance</tt>.
193
194  If \c n is 0 or 1, this function will call the error handler.
195  */
196  template<class vec_t> double vector_variance(const vec_t &data) {
197  return vector_variance(data.size(),data);
198  }
199
200  /** \brief Standard deviation with specified mean known in advance
201
202  This function computes
203  \f[
204  \sqrt{\frac{1}{N} \sum_{i} \left( x_i - \mu \right)^2}
205  \f]
206  where the value of \f$\mu \f$ is given in \c mean.
207
208  This function produces the same results
209  as <tt>gsl_stats_sd_with_fixed_mean()</tt>.
210
211  If \c n is zero, this function will return zero without calling
212  the error handler.
213  */
214  template<class vec_t>
215  double vector_stddev_fmean(size_t n, const vec_t &data,
216  double mean) {
217  double sd=vector_variance_fmean<vec_t>(n,data,mean);
218  return std::sqrt(sd);
219  }
220
221  /** \brief Standard deviation with specified mean known in advance
222
223  This function computes
224  \f[
225  \sqrt{\frac{1}{N} \sum_{i} \left( x_i - \mu \right)^2}
226  \f]
227  where the value of \f$\mu \f$ is given in \c mean.
228
229  This function produces the same results
230  as <tt>gsl_stats_sd_with_fixed_mean()</tt>.
231
232  If \c n is zero, this function will return zero without calling
233  the error handler.
234  */
235  template<class vec_t>
236  double vector_stddev_fmean(const vec_t &data, double mean) {
237  return vector_stddev_fmean(data.size(),data,mean);
238  }
239
240  /** \brief Standard deviation with specified mean
241
242  This function computes
243  \f[
244  \sqrt{\frac{1}{N-1} \sum_{i} \left( x_i - \mu \right)^2}
245  \f]
246  where \f$\mu \f$ is the mean computed with \ref vector_mean().
247
248  This function produces the same results
249  as <tt>gsl_stats_sd()</tt>.
250
251  If \c n is 0 or 1, this function will call the error handler.
252  */
253  template<class vec_t> double vector_stddev(size_t n, const vec_t &data) {
254
255  if (n<2) {
256  O2SCL_ERR2("Cannot compute std. dev. with less than 2 elements",
257  " in vector_stddev().",exc_einval);
258  }
259
260  double mean=vector_mean<vec_t>(n,data);
261  double var=vector_variance_fmean<vec_t>(n,data,mean);
262  return std::sqrt(var*n/(n-1));
263  }
264
265  /** \brief Standard deviation with specified mean
266
267  This function computes
268  \f[
269  \sqrt{\frac{1}{N-1} \sum_{i} \left( x_i - \mu \right)^2}
270  \f]
271  where \f$\mu \f$ is the mean computed with \ref vector_mean().
272
273  This function produces the same results
274  as <tt>gsl_stats_sd()</tt>.
275
276  If \c n is 0 or 1, this function will call the error handler.
277  */
278  template<class vec_t> double vector_stddev(const vec_t &data) {
279  return vector_stddev(data.size(),data);
280  }
281
282  /** \brief Standard deviation with specified mean
283
284  This function computes
285  \f[
286  \sqrt{\frac{1}{N-1} \sum_{i} \left( x_i - \mu \right)^2}
287  \f]
288  where the value of \f$\mu \f$ is given in \c mean.
289
290  This function produces the same results
291  as <tt>gsl_stats_sd_m()</tt>.
292
293  If \c n is 0 or 1, this function will call the error
294  handler.
295  */
296  template<class vec_t> double vector_stddev(size_t n, const vec_t &data,
297  double mean) {
298
299  if (n<2) {
300  O2SCL_ERR2("Cannot compute std. dev. with less than 2 elements",
301  " in vector_stddev().",exc_einval);
302  }
303
304  double sd=vector_variance_fmean<vec_t>(n,data,mean);
305  return std::sqrt(sd*n/(n-1));
306  }
307
308  /** \brief Standard deviation with specified mean
309
310  This function computes
311  \f[
312  \sqrt{\frac{1}{N-1} \sum_{i} \left( x_i - \mu \right)^2}
313  \f]
314  where the value of \f$\mu \f$ is given in \c mean.
315
316  This function produces the same results
317  as <tt>gsl_stats_sd_m()</tt>.
318
319  If \c n is 0 or 1, this function will call the error
320  handler.
321  */
322  template<class vec_t> double vector_stddev(const vec_t &data, double mean) {
323  return vector_stddev(data.size(),data,mean);
324  }
325
326  /** \brief Absolute deviation from the specified mean
327
328  This function computes
329  \f[
330  \sum_i | x_i - \mu |
331  \f]
332  where the value of \f$\mu \f$ is given in \c mean.
333
334  This function produces the same results
335  as <tt>gsl_stats_absdev_m()</tt>.
336
337  If \c n is zero, this function will return zero
338  without calling the error handler.
339  */
340  template<class vec_t> double vector_absdev(size_t n, const vec_t &data,
341  double mean) {
342
343  if (n==0) return 0.0;
344
345  long double sum=0.0;
346  for(size_t i=0;i<n;i++) {
347  sum+=fabs(data[i]-mean);
348  }
349  return sum/n;
350  }
351
352  /** \brief Absolute deviation from the specified mean
353
354  This function computes
355  \f[
356  \sum_i | x_i - \mu |
357  \f]
358  where the value of \f$\mu \f$ is given in \c mean.
359
360  This function produces the same results
361  as <tt>gsl_stats_absdev_m()</tt>.
362
363  If \c n is zero, this function will return zero
364  without calling the error handler.
365  */
366  template<class vec_t> double vector_absdev(const vec_t &data,
367  double mean) {
368  return vector_absdev(data.size(),data,mean);
369  }
370
371  /** \brief Absolute deviation from the computed mean
372
373  This function computes
374  \f[
375  \sum_i | x_i - \mu |
376  \f]
377  where the value of \f$\mu \f$ is mean as computed
378  from \ref vector_mean().
379
380  This function produces the same results
381  as <tt>gsl_stats_absdev()</tt>.
382
383  If \c n is zero, this function will return zero
384  without calling the error handler.
385  */
386  template<class vec_t>
387  double vector_absdev(size_t n, const vec_t &data) {
388  double mean=vector_mean<vec_t>(n,data);
389  return vector_absdev(n,data,mean);
390  }
391
392  /** \brief Absolute deviation from the computed mean
393
394  This function computes
395  \f[
396  \sum_i | x_i - \mu |
397  \f]
398  where the value of \f$\mu \f$ is mean as computed
399  from \ref vector_mean().
400
401  This function produces the same results
402  as <tt>gsl_stats_absdev()</tt>.
403
404  If \c n is zero, this function will return zero
405  without calling the error handler.
406  */
407  template<class vec_t>
408  double vector_absdev(const vec_t &data) {
409  return vector_absdev(data.size(),data);
410  }
411
412  /** \brief Skewness with specified mean and standard deviation
413
414  This function computes
415  \f[
416  \frac{1}{N} \sum_i \left[
417  \frac{ \left(x_i - \mu \right)}{ \sigma }\right]^3
418  \f]
419  where the values of \f$\mu \f$ and \f$\sigma \f$
420  are given in \c mean and \c stddev.
421
422  This function produces the same results
423  as <tt>gsl_stats_skew_m_sd()</tt>.
424
425  If \c n is zero, this function will return zero
426  without calling the error handler.
427  */
428  template<class vec_t> double vector_skew(size_t n, const vec_t &data,
429  double mean, double stddev) {
430  long double skew=0.0;
431  for(size_t i=0;i<n;i++) {
432  long double x=(data[i]-mean)/stddev;
433  skew+=(x*x*x-skew)/(i+1);
434  }
435  return skew;
436  }
437
438  /** \brief Skewness with specified mean and standard deviation
439
440  This function computes
441  \f[
442  \frac{1}{N} \sum_i \left[
443  \frac{ \left(x_i - \mu \right)}{ \sigma }\right]^3
444  \f]
445  where the values of \f$\mu \f$ and \f$\sigma \f$
446  are given in \c mean and \c stddev.
447
448  This function produces the same results
449  as <tt>gsl_stats_skew_m_sd()</tt>.
450
451  If \c n is zero, this function will return zero
452  without calling the error handler.
453  */
454  template<class vec_t> double vector_skew(const vec_t &data,
455  double mean, double stddev) {
456  return vector_skew(data.size(),data,mean,stddev);
457  }
458
459  /** \brief Skewness with computed mean and standard deviation
460
461  This function computes
462  \f[
463  \frac{1}{N} \sum_i \left[
464  \frac{ \left(x_i - \mu \right)}{ \sigma }\right]^3
465  \f]
466  where the values of \f$\mu \f$ and \f$\sigma \f$
467  are computed using \ref vector_mean() and \ref vector_stddev().
468
469  This function produces the same results
470  as <tt>gsl_stats_skew()</tt>.
471
472  If \c n is zero, this function will return zero
473  without calling the error handler.
474  */
475  template<class vec_t> double vector_skew(size_t n, const vec_t &data) {
476  double mean=vector_mean<vec_t>(n,data);
477  double sd=vector_stddev<vec_t>(n,data,mean);
478  return vector_skew(n,data,mean,sd);
479  }
480
481  /** \brief Skewness with computed mean and standard deviation
482
483  This function computes
484  \f[
485  \frac{1}{N} \sum_i \left[
486  \frac{ \left(x_i - \mu \right)}{ \sigma }\right]^3
487  \f]
488  where the values of \f$\mu \f$ and \f$\sigma \f$
489  are computed using \ref vector_mean() and \ref vector_stddev().
490
491  This function produces the same results
492  as <tt>gsl_stats_skew()</tt>.
493
494  If \c n is zero, this function will return zero
495  without calling the error handler.
496  */
497  template<class vec_t> double vector_skew(const vec_t &data) {
498  return vector_skew(data.size(),data);
499  }
500
501  /** \brief Kurtosis with specified mean and standard deviation
502
503  This function computes
504  \f[
505  -3 + \frac{1}{N} \sum_i \left[
506  \frac{ \left(x_i - \mu \right)}{ \sigma }\right]^4
507  \f]
508  where the values of \f$\mu \f$ and \f$\sigma \f$
509  are given in \c mean and \c stddev.
510
511  This function produces the same results
512  as <tt>gsl_stats_kurtosis_m_sd()</tt>.
513
514  If \c n is zero, this function will return zero
515  without calling the error handler.
516  */
517  template<class vec_t>
518  double vector_kurtosis(size_t n, const vec_t &data, double mean,
519  double stddev) {
520  long double avg=0.0;
521  for(size_t i=0;i<n;i++) {
522  long double x=(data[i]-mean)/stddev;
523  avg+=(x*x*x*x-avg)/(i+1);
524  }
525  return avg-3.0;
526  }
527
528  /** \brief Kurtosis with specified mean and standard deviation
529
530  This function computes
531  \f[
532  -3 + \frac{1}{N} \sum_i \left[
533  \frac{ \left(x_i - \mu \right)}{ \sigma }\right]^4
534  \f]
535  where the values of \f$\mu \f$ and \f$\sigma \f$
536  are given in \c mean and \c stddev.
537
538  This function produces the same results
539  as <tt>gsl_stats_kurtosis_m_sd()</tt>.
540
541  If \c n is zero, this function will return zero
542  without calling the error handler.
543  */
544  template<class vec_t>
545  double vector_kurtosis(const vec_t &data, double mean,
546  double stddev) {
547  return vector_kurtosis(data.size(),data,mean,stddev);
548  }
549
550  /** \brief Kurtosis with computed mean and standard deviation
551
552  This function computes
553  \f[
554  -3 + \frac{1}{N} \sum_i \left[
555  \frac{ \left(x_i - \mu \right)}{ \sigma }\right]^4
556  \f]
557  where the values of \f$\mu \f$ and \f$\sigma \f$
558  are computed using \ref vector_mean() and \ref vector_stddev().
559
560  This function produces the same results
561  as <tt>gsl_stats_kurtosis()</tt>.
562
563  If \c n is zero, this function will return zero
564  without calling the error handler.
565  */
566  template<class vec_t> double vector_kurtosis(size_t n, const vec_t &data) {
567  double mean=vector_mean<vec_t>(n,data);
568  double sd=vector_stddev<vec_t>(n,data,mean);
569  return vector_kurtosis(n,data,mean,sd);
570  }
571
572  /** \brief Kurtosis with computed mean and standard deviation
573
574  This function computes
575  \f[
576  -3 + \frac{1}{N} \sum_i \left[
577  \frac{ \left(x_i - \mu \right)}{ \sigma }\right]^4
578  \f]
579  where the values of \f$\mu \f$ and \f$\sigma \f$
580  are computed using \ref vector_mean() and \ref vector_stddev().
581
582  This function produces the same results
583  as <tt>gsl_stats_kurtosis()</tt>.
584
585  If \c n is zero, this function will return zero
586  without calling the error handler.
587  */
588  template<class vec_t> double vector_kurtosis(const vec_t &data) {
589  return vector_kurtosis(data.size(),data);
590  }
591
592  /** \brief Lag-1 autocorrelation
593
594  This function computes
595  \f[
596  \left[
597  \sum_i \left(x_i - \mu\right) \left(x_{i-1} - \mu \right)
598  \right] \left[
599  \sum_i \left(x_i - \mu\right)^2
600  \right]^{-1}
601  \f]
602
603  This function produces the same results
604  as <tt>gsl_stats_lag1_autocorrelation_m()</tt>.
605
606  If \c n is less than 2, this function will call the error handler.
607  */
608  template<class vec_t>
609  double vector_lag1_autocorr(size_t n, const vec_t &data, double mean) {
610
611  if (n<2) {
612  O2SCL_ERR2("Cannot compute lag1 with less than 2 elements",
613  " in vector_lag1_autocorr().",exc_einval);
614  }
615
616  long double q=0.0;
617  long double v=(data[0]-mean)*(data[0]-mean);
618  for(size_t i=1;i<n;i++) {
619  long double delta0=data[i-1]-mean;
620  long double delta1=data[i]-mean;
621  q+=(delta0*delta1-q)/(i+1);
622  v+=(delta1*delta1-v)/(i+1);
623  }
624
625  return q/v;
626  }
627
628  /** \brief Lag-1 autocorrelation
629
630  This function computes
631  \f[
632  \left[
633  \sum_i \left(x_i - \mu\right) \left(x_{i-1} - \mu \right)
634  \right] \left[
635  \sum_i \left(x_i - \mu\right)^2
636  \right]^{-1}
637  \f]
638
639  This function produces the same results
640  as <tt>gsl_stats_lag1_autocorrelation_m()</tt>.
641
642  If \c n is less than 2, this function will call the error handler.
643  */
644  template<class vec_t>
645  double vector_lag1_autocorr(const vec_t &data, double mean) {
646  return vector_lag1_autocorr(data.size(),data,mean);
647  }
648
649  /** \brief Lag-1 autocorrelation
650
651  This function computes
652  \f[
653  \left[
654  \sum_i \left(x_i - \mu\right) \left(x_{i-1} - \mu \right)
655  \right] \left[
656  \sum_i \left(x_i - \mu\right)^2
657  \right]^{-1}
658  \f]
659
660  This function produces the same results
661  as <tt>gsl_stats_lag1_autocorrelation()</tt>.
662
663  If \c n is less than 2, this function will call the error handler.
664  */
665  template<class vec_t> double vector_lag1_autocorr
666  (size_t n, const vec_t &data) {
667  double mean=vector_mean<vec_t>(n,data);
668  return vector_lag1_autocorr(n,data,mean);
669  }
670
671  /** \brief Lag-1 autocorrelation
672
673  This function computes
674  \f[
675  \left[
676  \sum_i \left(x_i - \mu\right) \left(x_{i-1} - \mu \right)
677  \right] \left[
678  \sum_i \left(x_i - \mu\right)^2
679  \right]^{-1}
680  \f]
681
682  This function produces the same results
683  as <tt>gsl_stats_lag1_autocorrelation()</tt>.
684
685  If \c n is less than 2, this function will call the error handler.
686  */
687  template<class vec_t> double vector_lag1_autocorr(const vec_t &data) {
688  return vector_lag1_autocorr(data.size(),data);
689  }
690
691  /** \brief Lag-k autocorrelation
692
693  This function computes
694  \f[
695  \left[
696  \sum_i \left(x_i - \mu\right) \left(x_{i-k} - \mu \right)
697  \right] \left[
698  \sum_i \left(x_i - \mu\right)^2
699  \right]^{-1}
700  \f]
701
702  If <tt>n<=k</tt>, this function will call the error handler.
703  */
704  template<class vec_t>
705  double vector_lagk_autocorr(size_t n, const vec_t &data, size_t k,
706  double mean) {
707
708  if (n<=k) {
709  O2SCL_ERR2("Not enough elements ",
710  " in vector_lagk_autocorr().",exc_einval);
711  }
712
713  long double q=0.0, v=0.0;
714  for(size_t i=0;i<k;i++) {
715  q+=0.0;
716  v+=(data[i]-mean)*(data[i]-mean)/(i+1);
717  }
718  for(size_t i=k;i<n;i++) {
719  long double delta0=data[i-k]-mean;
720  long double delta1=data[i]-mean;
721  q+=(delta0*delta1-q)/(i+1);
722  v+=(delta1*delta1-v)/(i+1);
723  }
724  return q/v;
725  }
726
727  /** \brief Lag-k autocorrelation
728
729  This function computes
730  \f[
731  \left[
732  \sum_i \left(x_i - \mu\right) \left(x_{i-k} - \mu \right)
733  \right] \left[
734  \sum_i \left(x_i - \mu\right)^2
735  \right]^{-1}
736  \f]
737
738  If <tt>n<=k</tt>, this function will call the error handler.
739  */
740  template<class vec_t>
741  double vector_lagk_autocorr(const vec_t &data, size_t k,
742  double mean) {
743  return vector_lagk_autocorr(data.size(),k,mean);
744  }
745
746  /** \brief Lag-k autocorrelation
747
748  This function computes
749  \f[
750  \left[
751  \sum_i \left(x_i - \mu\right) \left(x_{i-k} - \mu \right)
752  \right] \left[
753  \sum_i \left(x_i - \mu\right)^2
754  \right]^{-1}
755  \f]
756
757  If <tt>n<=k</tt>, this function will call the error handler.
758  */
759  template<class vec_t> double vector_lagk_autocorr
760  (size_t n, const vec_t &data, size_t k) {
761  double mean=vector_mean<vec_t>(n,data);
762  return vector_lagk_autocorr(n,data,k,mean);
763  }
764
765  /** \brief Lag-k autocorrelation
766
767  This function computes
768  \f[
769  \left[
770  \sum_i \left(x_i - \mu\right) \left(x_{i-k} - \mu \right)
771  \right] \left[
772  \sum_i \left(x_i - \mu\right)^2
773  \right]^{-1}
774  \f]
775
776  If <tt>n<=k</tt>, this function will call the error handler.
777  */
778  template<class vec_t> double vector_lagk_autocorr
779  (const vec_t &data, size_t k) {
780  return vector_lagk_autocorr(data.size(),data,k);
781  }
782
783  /** \brief Compute the covariance of two vectors
784
785  This function computes
786  \f[
787  \frac{1}{n-1} \sum_i \left(x_i - {\bar{x}}\right)
788  \left(y_i - {\bar{y}}\right)
789  \f]
790  where \f${\bar{x}} \f$ and \f${\bar{y}} \f$ are specified
791  in \c mean1 and \c mean2, respectively.
792
793  This function produces the same results
794  as <tt>gsl_stats_covariance_m()</tt>.
795
796  If \c n is zero, this function will return zero
797  without calling the error handler.
798  */
799  template<class vec_t, class vec2_t>
800  double vector_covariance(size_t n, const vec_t &data1, const vec2_t &data2,
801  double mean1, double mean2) {
802  double covar=0.0;
803  for(size_t i=0;i<n;i++) {
804  double delta1=(data1[i]-mean1);
805  double delta2=(data2[i]-mean2);
806  covar+=(delta1*delta2-covar)/(i+1);
807  }
808  return covar*n/(n-1);
809  }
810
811  /** \brief Compute the covariance of two vectors
812
813  This function computes
814  \f[
815  \frac{1}{n-1} \sum_i \left(x_i - {\bar{x}}\right)
816  \left(y_i - {\bar{y}}\right)
817  \f]
818  where \f${\bar{x}} \f$ and \f${\bar{y}} \f$ are specified
819  in \c mean1 and \c mean2, respectively.
820
821  This function produces the same results
822  as <tt>gsl_stats_covariance_m()</tt>.
823
824  If \c n is zero, this function will return zero
825  without calling the error handler.
826  */
827  template<class vec_t, class vec2_t>
828  double vector_covariance(const vec_t &data1, const vec2_t &data2,
829  double mean1, double mean2) {
830  return vector_covariance(data1.size(),data1,data2,mean1,mean2);
831  }
832
833  /** \brief Compute the covariance of two vectors
834
835  This function computes
836  \f[
837  \frac{1}{n-1} \sum_i \left(x_i - {\bar{x}}\right)
838  \left(y_i - {\bar{y}}\right)
839  \f]
840  where \f${\bar{x}} \f$ and \f${\bar{y}} \f$ are
841  the averages of \c data1 and \c data2 and are computed
842  automatically using \ref vector_mean().
843
844  This function produces the same
845  results as <tt>gsl_stats_covariance()</tt>.
846
847  If \c n is zero, this function will return zero
848  without calling the error handler.
849  */
850  template<class vec_t, class vec2_t>
851  double vector_covariance(size_t n, const vec_t &data1,
852  const vec2_t &data2) {
853  double covar=0.0;
854  double mean1=vector_mean<vec_t>(n,data1);
855  double mean2=vector_mean<vec_t>(n,data2);
856  for(size_t i=0;i<n;i++) {
857  long double delta1=(data1[i]-mean1);
858  long double delta2=(data2[i]-mean2);
859  covar+=(delta1*delta2-covar)/(i+1);
860  }
861  return covar*n/(n-1);
862  }
863
864  /** \brief Compute the covariance of two vectors
865
866  This function computes
867  \f[
868  \frac{1}{n-1} \sum_i \left(x_i - {\bar{x}}\right)
869  \left(y_i - {\bar{y}}\right)
870  \f]
871  where \f${\bar{x}} \f$ and \f${\bar{y}} \f$ are
872  the averages of \c data1 and \c data2 and are computed
873  automatically using \ref vector_mean().
874
875  This function produces the same
876  results as <tt>gsl_stats_covariance()</tt>.
877
878  If \c n is zero, this function will return zero
879  without calling the error handler.
880  */
881  template<class vec_t, class vec2_t>
882  double vector_covariance(const vec_t &data1,
883  const vec2_t &data2) {
884  return vector_covariance(data1.size(),data1,data2);
885  }
886
887  /** \brief Pearson's correlation
888
889  This function computes the Pearson correlation coefficient
890  between \c data1 and \c data2 .
891
892  This function produces the same
893  results as <tt>gsl_stats_correlation()</tt>.
894
895  \comment
896  r = cov(x, y) / (\Hat\sigma_x \Hat\sigma_y)
897  = {1/(n-1) \sum (x_i - \Hat x) (y_i - \Hat y)
898  \over
899  \sqrt{1/(n-1) \sum (x_i - \Hat x)^2} \sqrt{1/(n-1)
900  \sum (y_i - \Hat y)^2}
901  }
902  \endcomment
903
904  If \c n is zero, this function will call the error handler.
905  */
906  template<class vec_t, class vec2_t>
907  double vector_correlation(size_t n, const vec_t &data1,
908  const vec2_t &data2) {
909  size_t i;
910
911  if (n<1) {
912  O2SCL_ERR2("Cannot compute correlation with no elements",
913  " in vector_correlation().",exc_einval);
914  }
915
916  double sum_xsq=0.0;
917  double sum_ysq=0.0;
918  double sum_cross=0.0;
919  double ratio;
920  double delta_x, delta_y;
921  double mean_x, mean_y;
922  double r;
923
924  /*
925  * Compute:
926  * sum_xsq = Sum [ (x_i - mu_x)^2 ],
927  * sum_ysq = Sum [ (y_i - mu_y)^2 ] and
928  * sum_cross = Sum [ (x_i - mu_x) * (y_i - mu_y) ]
929  * using the above relation from Welford's paper
930  */
931
932  mean_x=data1[0];
933  mean_y=data2[0];
934
935  for (i=1; i < n; ++i) {
936  ratio=i / (i + 1.0);
937  delta_x=data1[i] - mean_x;
938  delta_y=data2[i] - mean_y;
939  sum_xsq += delta_x * delta_x * ratio;
940  sum_ysq += delta_y * delta_y * ratio;
941  sum_cross += delta_x * delta_y * ratio;
942  mean_x += delta_x / (i + 1.0);
943  mean_y += delta_y / (i + 1.0);
944  }
945
946  r=sum_cross / (std::sqrt(sum_xsq) * std::sqrt(sum_ysq));
947
948  return r;
949  }
950
951  /** \brief Pearson's correlation
952
953  This function computes the Pearson correlation coefficient
954  between \c data1 and \c data2 .
955
956  This function produces the same
957  results as <tt>gsl_stats_correlation()</tt>.
958
959  \comment
960  r = cov(x, y) / (\Hat\sigma_x \Hat\sigma_y)
961  = {1/(n-1) \sum (x_i - \Hat x) (y_i - \Hat y)
962  \over
963  \sqrt{1/(n-1) \sum (x_i - \Hat x)^2} \sqrt{1/(n-1)
964  \sum (y_i - \Hat y)^2}
965  }
966  \endcomment
967
968  If \c n is zero, this function will call the error handler.
969  */
970  template<class vec_t, class vec2_t>
971  double vector_correlation(const vec_t &data1,
972  const vec2_t &data2) {
973  return vector_correlation(data1.size(),data1,data2);
974  }
975
976  /** \brief The pooled variance of two vectors
977
978  This function computes
979  \f[
980  s_{p}^2 = \frac{(n_1-1)s_1^2+(n_2-1)s_2^2}{n_1+n_2-2}
981  \f]
982  where \f$n_i \f$ is the number of elements in vector \f$i \f$
983  and \f$s_i^2 \f$ is the variance of vector \f$i \f$.
984
985  From http://en.wikipedia.org/wiki/Pooled_variance, "Under the
986  assumption of equal population variances, the pooled sample
987  variance provides a higher precision estimate of variance than
988  the individual sample variances."
989
990  This function produces the same
991  results as <tt>gsl_stats_pvariance()</tt>.
992  */
993  template<class vec_t, class vec2_t>
994  double vector_pvariance(size_t n1, const vec_t &data1,
995  size_t n2, const vec2_t &data2) {
996  double var1=vector_variance<vec_t>(n1,data1);
997  double var2=vector_variance<vec2_t>(n2,data2);
998  return (((n1-1)*var1)+((n2-1)*var2))/(n1+n2-2);
999  }
1000
1001  /** \brief The pooled variance of two vectors
1002
1003  This function computes
1004  \f[
1005  s_{p}^2 = \frac{(n_1-1)s_1^2+(n_2-1)s_2^2}{n_1+n_2-2}
1006  \f]
1007  where \f$n_i \f$ is the number of elements in vector \f$i \f$
1008  and \f$s_i^2 \f$ is the variance of vector \f$i \f$.
1009
1010  From http://en.wikipedia.org/wiki/Pooled_variance, "Under the
1011  assumption of equal population variances, the pooled sample
1012  variance provides a higher precision estimate of variance than
1013  the individual sample variances."
1014
1015  This function produces the same
1016  results as <tt>gsl_stats_pvariance()</tt>.
1017  */
1018  template<class vec_t, class vec2_t>
1019  double vector_pvariance(const vec_t &data1,
1020  const vec2_t &data2) {
1021  return vector_pvariance(data1.size(),data1,data2.size(),data2);
1022  }
1023
1024  /** \brief Quantile from sorted data (ascending only)
1025
1026  This function returns the quantile \c f of data which
1027  has already been sorted in ascending order. The quantile,
1028  \f$q \f$ , is
1029  found by interpolation using
1030  \f[
1031  q = \left(1-\delta\right) x_i \delta x_{i+1}
1032  \f]
1033  where \f$i = \mathrm{floor}[ (n-1)f ] \f$ and
1034  \f$\delta = (n-1)f -i \f$ .
1035
1036  This function produces the same
1037  results as <tt>gsl_stats_quantile_from_sorted_data()</tt>.
1038
1039  No checks are made to ensure the data is sorted, or to ensure
1040  that \f$0 \leq 0 \leq 1 \f$. If \c n is zero, this function
1041  will return zero without calling the error handler.
1042  */
1043  template<class vec_t>
1044  double vector_quantile_sorted(size_t n, const vec_t &data,
1045  const double f) {
1046
1047  double index=f*(n-1);
1048  size_t lhs=((size_t)index);
1049  double delta=index-lhs;
1050  if (n==0) return 0.0;
1051  if (lhs==n-1) return data[lhs];
1052  return (1-delta)*data[lhs]+delta*data[lhs+1];
1053  }
1054
1055  /** \brief Quantile from sorted data (ascending only)
1056
1057  This function returns the quantile \c f of data which
1058  has already been sorted in ascending order. The quantile,
1059  \f$q \f$ , is
1060  found by interpolation using
1061  \f[
1062  q = \left(1-\delta\right) x_i \delta x_{i+1}
1063  \f]
1064  where \f$i = \mathrm{floor}[ (n-1)f ] \f$ and
1065  \f$\delta = (n-1)f -i \f$ .
1066
1067  This function produces the same
1068  results as <tt>gsl_stats_quantile_from_sorted_data()</tt>.
1069
1070  No checks are made to ensure the data is sorted, or to ensure
1071  that \f$0 \leq 0 \leq 1 \f$. If \c n is zero, this function
1072  will return zero without calling the error handler.
1073  */
1074  template<class vec_t>
1075  double vector_quantile_sorted(const vec_t &data, const double f) {
1076  return vector_quantile_sorted<vec_t>(data.size(),data,f);
1077  }
1078
1079  /** \brief Return the median of sorted (ascending or descending) data
1080
1081  This function returns the median of sorted data (either
1082  ascending or descending), assuming the data has already been
1083  sorted. When the data set has an odd number of elements, the
1084  median is the value of the element at index \f$(n-1)/2 \f$,
1085  otherwise, the median is taken to be the average of the elements
1086  at indices \f$(n-1)/2 \f$ and \f$n/2 \f$ .
1087
1088  This function produces the same
1089  results as <tt>gsl_stats_median_from_sorted_data()</tt>.
1090
1091  No checks are made to ensure the data is sorted. If \c n is
1092  zero, this function will return zero without calling the error
1093  handler.
1094  */
1095  template<class vec_t>
1096  double vector_median_sorted(size_t n, const vec_t &data) {
1097
1098  if (n==0) return 0.0;
1099
1100  size_t lhs=(n-1)/2;
1101  size_t rhs=n/2;
1102
1103  if (lhs==rhs) return data[lhs];
1104
1105  return (data[lhs]+data[rhs])/2.0;
1106  }
1107
1108  /** \brief Return the median of sorted (ascending or descending) data
1109
1110  This function returns the median of sorted data (either
1111  ascending or descending), assuming the data has already been
1112  sorted. When the data set has an odd number of elements, the
1113  median is the value of the element at index \f$(n-1)/2 \f$,
1114  otherwise, the median is taken to be the average of the elements
1115  at indices \f$(n-1)/2 \f$ and \f$n/2 \f$ .
1116
1117  This function produces the same
1118  results as <tt>gsl_stats_median_from_sorted_data()</tt>.
1119
1120  No checks are made to ensure the data is sorted. If \c n is
1121  zero, this function will return zero without calling the error
1122  handler.
1123  */
1124  template<class vec_t>
1125  double vector_median_sorted(const vec_t &data) {
1126  return vector_median_sorted<vec_t>(data.size(),data);
1127  }
1128
1129  /** \brief Compute the chi-squared statistic
1130
1131  This function computes
1132  \f[
1133  \sum_i \left( \frac{\mathrm{obs}_i - \mathrm{exp}_i}
1134  {\mathrm{err}_i}\right)^2
1135  \f]
1136  where \f$\mathrm{obs} \f$ are the observed values,
1137  \f$\mathrm{exp} \f$ are the expected values, and
1138  \f$\mathrm{err} \f$ are the errors.
1139  */
1140  template<class vec_t, class vec2_t, class vec3_t>
1141  double vector_chi_squared(size_t n, const vec_t &obs, const vec2_t &exp,
1142  const vec3_t &err) {
1143  double chi2=0.0;
1144  for(size_t i=0;i<n;i++) {
1145  chi2+=pow((obs[i]-exp[i])/err[i],2.0);
1146  }
1147  return chi2;
1148  }
1149
1150  /** \brief Compute the chi-squared statistic
1151
1152  This function computes
1153  \f[
1154  \sum_i \left( \frac{\mathrm{obs}_i - \mathrm{exp}_i}
1155  {\mathrm{err}_i}\right)^2
1156  \f]
1157  where \f$\mathrm{obs} \f$ are the observed values,
1158  \f$\mathrm{exp} \f$ are the expected values, and
1159  \f$\mathrm{err} \f$ are the errors.
1160  */
1161  template<class vec_t, class vec2_t, class vec3_t>
1162  double vector_chi_squared(const vec_t &obs, const vec2_t &exp,
1163  const vec3_t &err) {
1164  return vector_chi_squared<vec_t,vec2_t,vec3_t>(obs.size(),obs,exp,err);
1165  }
1166  //@}
1167
1168  /// \name Weighted vector functions
1169  //@{
1170  /** \brief Compute the mean of weighted data
1171
1172  This function computes
1173  \f[
1174  \left( \sum_i w_i x_i \right) \left( \sum_i w_i \right)^{-1}
1175  \f]
1176
1177  This function produces the same results
1178  as <tt>gsl_stats_wmean()</tt>.
1179
1180  \comment
1181  M(n) = M(n-1) + (data[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
1182  W(n) = W(n-1) + w(n)
1183  \endcomment
1184  */
1185  template<class vec_t, class vec2_t>
1186  double wvector_mean(size_t n, const vec_t &data, const vec2_t &weights) {
1187
1188  long double wmean=0.0;
1189  long double W=0.0;
1190  for(size_t i=0;i<n;i++) {
1191  double wi=weights[i];
1192  if (wi>0.0) {
1193  W+=wi;
1194  wmean+=(data[i]-wmean)*(wi/W);
1195  }
1196  }
1197
1198  return wmean;
1199  }
1200
1201  /** \brief Compute the mean of weighted data
1202
1203  This function computes
1204  \f[
1205  \left( \sum_i w_i x_i \right) \left( \sum_i w_i \right)^{-1}
1206  \f]
1207
1208  This function produces the same results
1209  as <tt>gsl_stats_wmean()</tt>.
1210
1211  \comment
1212  M(n) = M(n-1) + (data[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
1213  W(n) = W(n-1) + w(n)
1214  \endcomment
1215  */
1216  template<class vec_t, class vec2_t>
1217  double wvector_mean(const vec_t &data, const vec2_t &weights) {
1218  return wvector_mean<vec_t,vec2_t>(data.size(),data,weights);
1219  }
1220
1221  /** \brief Compute a normalization factor for weighted data
1222
1223  This function is used internally in \ref wvector_variance(size_t
1224  n, vec_t &data, const vec2_t &weights, double wmean) and \ref
1225  wvector_stddev(size_t n, vec_t &data, const vec2_t &weights, double
1226  wmean) .
1227  */
1228  template<class vec_t> double wvector_factor(size_t n, const vec_t &weights) {
1229
1230  long double a=0.0;
1231  long double b=0.0;
1232  long double factor;
1233  for(size_t i=0;i<n;i++) {
1234  double wi=weights[i];
1235  if (wi>0.0) {
1236  a+=wi;
1237  b+=wi*wi;
1238  }
1239  }
1240  factor=a*a/(a*a-b);
1241  return factor;
1242  }
1243
1244  /** \brief Compute a normalization factor for weighted data
1245
1246  This function is used internally in \ref wvector_variance(size_t
1247  n, vec_t &data, const vec2_t &weights, double wmean) and \ref
1248  wvector_stddev(size_t n, vec_t &data, const vec2_t &weights, double
1249  wmean) .
1250  */
1251  template<class vec_t> double wvector_factor(const vec_t &weights) {
1252  return wvector_factor<vec_t>(weights.size(),weights);
1253  }
1254
1255  /** \brief Compute the variance of a weighted vector with a mean
1257
1258  This function computes
1259  \f[
1260  \left[ \sum_i w_i \left(x_i-\mu\right)^2 \right]
1261  \left[ \sum_i w_i \right]^{-1}
1262  \f]
1263
1264  This function produces the same results
1265  as <tt>gsl_stats_wvariance_with_fixed_mean()</tt>.
1266
1267  */
1268  template<class vec_t, class vec2_t>
1269  double wvector_variance_fmean(size_t n, const vec_t &data,
1270  const vec2_t &weights, double wmean) {
1271  long double wvariance=0.0;
1272  long double W=0.0;
1273  for(size_t i=0;i<n;i++) {
1274  double wi=weights[i];
1275  if (wi>0.0) {
1276  const long double delta=data[i]-wmean;
1277  W+=wi;
1278  wvariance+=(delta*delta-wvariance)*(wi/W);
1279  }
1280  }
1281
1282  return wvariance;
1283  }
1284
1285  /** \brief Compute the variance of a weighted vector with a mean
1287
1288  This function computes
1289  \f[
1290  \left[ \sum_i w_i \left(x_i-\mu\right)^2 \right]
1291  \left[ \sum_i w_i \right]^{-1}
1292  \f]
1293
1294  This function produces the same results
1295  as <tt>gsl_stats_wvariance_with_fixed_mean()</tt>.
1296
1297  */
1298  template<class vec_t, class vec2_t>
1299  double wvector_variance_fmean(const vec_t &data,
1300  const vec2_t &weights, double wmean) {
1301  return wvector_variance_fmean(data.size(),data,weights,wmean);
1302  }
1303
1304  /** \brief Compute the variance of a weighted vector with
1305  specified mean
1306
1307  This function produces the same results
1308  as <tt>gsl_stats_wvariance_m()</tt>.
1309  */
1310  template<class vec_t, class vec2_t>
1311  double wvector_variance(size_t n, const vec_t &data,
1312  const vec2_t &weights, double wmean) {
1313
1314  const double variance=wvector_variance_fmean
1315  (n,data,weights,wmean);
1316  const double scale=wvector_factor(n,weights);
1317  const double wvar=scale*variance;
1318  return wvar;
1319  }
1320
1321  /** \brief Compute the variance of a weighted vector with
1322  specified mean
1323
1324  This function produces the same results
1325  as <tt>gsl_stats_wvariance_m()</tt>.
1326  */
1327  template<class vec_t, class vec2_t>
1328  double wvector_variance(const vec_t &data,
1329  const vec2_t &weights, double wmean) {
1330  return wvector_variance<vec_t,vec2_t>(data.size(),data,weights,wmean);
1331  }
1332
1333  /** \brief Compute the variance of a weighted vector where mean
1334  is computed automatically
1335
1336  This function produces the same results
1337  as <tt>gsl_stats_wvariance()</tt>.
1338  */
1339  template<class vec_t, class vec2_t>
1340  double wvector_variance(size_t n, const vec_t &data,
1341  const vec2_t &weights) {
1342
1343  double wmean=wvector_mean(n,data,weights);
1344  return wvector_variance<vec_t,vec2_t>(n,data,weights,wmean);
1345  }
1346
1347  /** \brief Compute the variance of a weighted vector where mean
1348  is computed automatically
1349
1350  This function produces the same results
1351  as <tt>gsl_stats_wvariance()</tt>.
1352  */
1353  template<class vec_t, class vec2_t>
1354  double wvector_variance(const vec_t &data, const vec2_t &weights) {
1355  return wvector_variance(data.size(),data,weights);
1356  }
1357
1358  /** \brief The weighted covariance of two vectors
1359
1360  \note Experimental
1361  */
1362  template<class vec_t, class vec2_t, class vec3_t>
1363  double wvector_covariance(size_t n, const vec_t &data1, const vec2_t &data2,
1364  const vec3_t &weights) {
1365  double mean1=wvector_mean(n,data1,weights);
1366  double mean2=wvector_mean(n,data2,weights);
1367  double covar=0.0;
1368  double W=0.0;
1369  for(size_t i=0;i<n;i++) {
1370  double wi=weights[i];
1371  if (wi>0.0) {
1372  W+=wi;
1373  double delta1=(data1[i]-mean1);
1374  double delta2=(data2[i]-mean2);
1375  covar+=(wi/W)*(delta1*delta2-covar);
1376  }
1377  }
1378  double scale=wvector_factor(n,weights);
1379  return covar*scale;
1380  }
1381
1382  /** \brief The weighted covariance of two vectors
1383
1384  \note Experimental
1385  */
1386  template<class vec_t, class vec2_t, class vec3_t>
1387  double wvector_covariance(const vec_t &data1, const vec2_t &data2,
1388  const vec3_t &weights) {
1389  return wvector_covariance<vec_t,vec2_t,vec3_t>
1390  (data1.size(),data1,data2,weights);
1391  }
1392
1393  /** \brief Compute the standard deviation of a weighted vector
1394  with a mean known in advance
1395
1396  This function produces the same results
1397  as <tt>gsl_stats_wsd_with_fixed_mean()</tt>.
1398  */
1399  template<class vec_t, class vec2_t>
1400  double wvector_stddev_fmean(size_t n, const vec_t &data,
1401  const vec2_t &weights, double wmean) {
1402  return sqrt(wvector_variance_fmean(n,data,weights,wmean));
1403  }
1404
1405  /** \brief Compute the standard deviation of a weighted vector
1406  with a mean known in advance
1407
1408  This function produces the same results
1409  as <tt>gsl_stats_wsd_with_fixed_mean()</tt>.
1410  */
1411  template<class vec_t, class vec2_t>
1412  double wvector_stddev_fmean(const vec_t &data,
1413  const vec2_t &weights, double wmean) {
1414  return wvector_stddev_fmean<vec_t,vec2_t>
1415  (data.size(),data,weights,wmean);
1416  }
1417
1418  /** \brief Compute the standard deviation of a weighted vector where mean
1419  is computed automatically
1420
1421  This function produces the same results
1422  as <tt>gsl_stats_wsd()</tt>.
1423  */
1424  template<class vec_t, class vec2_t>
1425  double wvector_stddev(size_t n, const vec_t &data,
1426  const vec2_t &weights) {
1427  double wmean=wvector_mean(n,data,weights);
1428  return sqrt(wvector_variance(n,data,weights,wmean));
1429  }
1430
1431  /** \brief Compute the standard deviation of a weighted vector where mean
1432  is computed automatically
1433
1434  This function produces the same results
1435  as <tt>gsl_stats_wsd()</tt>.
1436  */
1437  template<class vec_t, class vec2_t>
1438  double wvector_stddev(const vec_t &data, const vec2_t &weights) {
1439  return wvector_stddev(data.size(),data,weights);
1440  }
1441
1442  /** \brief Compute the standard deviation of a weighted vector with
1443  specified mean
1444
1445  This function produces the same results
1446  as <tt>gsl_stats_wsd_m()</tt>.
1447  */
1448  template<class vec_t, class vec2_t>
1449  double wvector_stddev(size_t n, const vec_t &data,
1450  const vec2_t &weights, double wmean) {
1451  const double variance=wvector_variance_fmean
1452  (n,data,weights,wmean);
1453  const double scale=wvector_factor(n,weights);
1454  double wvar=scale*variance;
1455  return sqrt(wvar);
1456  }
1457
1458  /** \brief Compute the standard deviation of a weighted vector with
1459  specified mean
1460
1461  This function produces the same results
1462  as <tt>gsl_stats_wsd_m()</tt>.
1463  */
1464  template<class vec_t, class vec2_t>
1465  double wvector_stddev(const vec_t &data,
1466  const vec2_t &weights, double wmean) {
1467  return wvector_stddev<vec_t,vec2_t>(data.size(),data,weights,wmean);
1468  }
1469
1470  /** \brief Compute the weighted sum of squares of data about the
1471  specified weighted mean
1472
1473  This function produces the same results
1474  as <tt>gsl_stats_wtss_m()</tt>.
1475  */
1476  template<class vec_t, class vec2_t>
1477  double wvector_sumsq(size_t n, const vec_t &data,
1478  const vec2_t &weights, double wmean) {
1479  long double wtss=0.0;
1480  for(size_t i=0;i<n;i++) {
1481  double wi=weights[i];
1482  if (wi>0.0) {
1483  const long double delta=data[i]-wmean;
1484  wtss+=wi*delta*delta;
1485  }
1486  }
1487
1488  return wtss;
1489  }
1490
1491  /** \brief Compute the weighted sum of squares of data about the
1492  specified weighted mean
1493
1494  This function produces the same results
1495  as <tt>gsl_stats_wtss_m()</tt>.
1496  */
1497  template<class vec_t, class vec2_t>
1498  double wvector_sumsq(const vec_t &data,
1499  const vec2_t &weights, double wmean) {
1500  return wvector_sumsq<vec_t,vec2_t>(data.size(),data,weights,wmean);
1501  }
1502
1503  /** \brief Compute the weighted sum of squares of data about the
1504  weighted mean
1505
1506  This function produces the same results
1507  as <tt>gsl_stats_wtss()</tt>.
1508  */
1509  template<class vec_t, class vec2_t>
1510  double wvector_sumsq(size_t n, const vec_t &data,
1511  const vec2_t &weights) {
1512
1513  double wmean=wvector_mean(n,data,weights);
1514  return wvector_sumsq(n,data,weights,wmean);
1515  }
1516
1517  /** \brief Compute the weighted sum of squares of data about the
1518  weighted mean
1519
1520  This function produces the same results
1521  as <tt>gsl_stats_wtss()</tt>.
1522  */
1523  template<class vec_t, class vec2_t>
1524  double wvector_sumsq(const vec_t &data, const vec2_t &weights) {
1525  return wvector_sumsq<vec_t,vec2_t>(data.size(),data,weights);
1526  }
1527
1528  /** \brief Compute the absolute deviation of data about a specified mean
1529
1530  This function produces the same results
1531  as <tt>gsl_stats_wabsdev_m()</tt>.
1532  */
1533  template<class vec_t, class vec2_t>
1534  double wvector_absdev(size_t n, const vec_t &data, const vec2_t &weights,
1535  double wmean) {
1536  long double wabsdev=0.0;
1537  long double W=0.0;
1538  for(size_t i=0;i<n;i++) {
1539  double wi=weights[i];
1540  if (wi>0.0) {
1541  const long double delta=fabs(data[i]-wmean);
1542  W+=wi;
1543  wabsdev+=(delta-wabsdev)*(wi/W);
1544  }
1545  }
1546  return wabsdev;
1547  }
1548
1549  /** \brief Compute the absolute deviation of data about a specified mean
1550
1551  This function produces the same results
1552  as <tt>gsl_stats_wabsdev_m()</tt>.
1553  */
1554  template<class vec_t, class vec2_t>
1555  double wvector_absdev(const vec_t &data, const vec2_t &weights,
1556  double wmean) {
1557  return wvector_absdev<vec_t,vec2_t>(data.size(),data,weights,wmean);
1558  }
1559
1560  /** \brief Compute the absolute deviation of data about a specified mean
1561
1562  This function produces the same results
1563  as <tt>gsl_stats_wabsdev()</tt>.
1564  */
1565  template<class vec_t, class vec2_t>
1566  double wvector_absdev(size_t n, const vec_t &data,
1567  const vec2_t &weights) {
1568
1569  double wmean=wvector_mean(n,data,weights);
1570  return wvector_absdev(n,data,weights,wmean);
1571  }
1572
1573  /** \brief Compute the absolute deviation of data about a specified mean
1574
1575  This function produces the same results
1576  as <tt>gsl_stats_wabsdev()</tt>.
1577  */
1578  template<class vec_t, class vec2_t>
1579  double wvector_absdev(const vec_t &data, const vec2_t &weights) {
1580  return wvector_absdev<vec_t,vec2_t>(data.size(),data,weights);
1581  }
1582
1583  /** \brief Compute the skewness of data with specified mean
1584  and standard deviation
1585
1586  This function produces the same results
1587  as <tt>gsl_stats_wskew_m_sd()</tt>.
1588  */
1589  template<class vec_t, class vec2_t>
1590  double wvector_skew(size_t n, const vec_t &data, const vec2_t &weights,
1591  double wmean, double wsd) {
1592  long double wskew=0.0;
1593  long double W=0.0;
1594  for(size_t i=0;i<n;i++) {
1595  double wi=weights[i];
1596  if (wi>0.0) {
1597  const long double x=(data[i]-wmean)/wsd;
1598  W+=wi;
1599  wskew+=(x*x*x-wskew)*(wi/W);
1600  }
1601  }
1602  return wskew;
1603  }
1604
1605  /** \brief Compute the skewness of data with specified mean
1606  and standard deviation
1607
1608  This function produces the same results
1609  as <tt>gsl_stats_wskew_m_sd()</tt>.
1610  */
1611  template<class vec_t, class vec2_t>
1612  double wvector_skew(const vec_t &data, const vec2_t &weights,
1613  double wmean, double wsd) {
1614  return wvector_skew<vec_t,vec2_t>(data.size(),data,weights,wmean,wsd);
1615  }
1616
1617  /** \brief Compute the skewness of data with specified mean
1618  and standard deviation
1619
1620  This function produces the same results
1621  as <tt>gsl_stats_wskew()</tt>.
1622  */
1623  template<class vec_t, class vec2_t>
1624  double wvector_skew(size_t n, const vec_t &data, const vec2_t &weights) {
1625  double wmean=wvector_mean(n,data,weights);
1626  double wsd=wvector_stddev(n,data,weights,wmean);
1627  return wvector_skew(n,data,weights,wmean,wsd);
1628  }
1629
1630  /** \brief Compute the skewness of data with specified mean
1631  and standard deviation
1632
1633  This function produces the same results
1634  as <tt>gsl_stats_wskew()</tt>.
1635  */
1636  template<class vec_t, class vec2_t>
1637  double wvector_skew(const vec_t &data, const vec2_t &weights) {
1638  return wvector_skew<vec_t,vec2_t>(data.size(),data,weights);
1639  }
1640
1641  /** \brief Compute the kurtosis of data with specified mean
1642  and standard deviation
1643
1644  This function produces the same results
1645  as <tt>gsl_stats_wkurtosis_m_sd()</tt>.
1646  */
1647  template<class vec_t, class vec2_t>
1648  double wvector_kurtosis(size_t n, const vec_t &data, const vec2_t &weights,
1649  double wmean, double wsd) {
1650  long double wavg=0.0;
1651  long double W=0.0;
1652  for(size_t i=0;i<n;i++) {
1653  double wi=weights[i];
1654  if (wi>0.0) {
1655  const long double x=(data[i]-wmean)/wsd;
1656  W+=wi;
1657  wavg+=(x*x*x*x-wavg)*(wi/W);
1658  }
1659  }
1660  return wavg-3.0;
1661  }
1662
1663  /** \brief Compute the kurtosis of data with specified mean
1664  and standard deviation
1665
1666  This function produces the same results
1667  as <tt>gsl_stats_wkurtosis_m_sd()</tt>.
1668  */
1669  template<class vec_t, class vec2_t>
1670  double wvector_kurtosis(const vec_t &data, const vec2_t &weights,
1671  double wmean, double wsd) {
1672  return wvector_kurtosis<vec_t,vec2_t>
1673  (data.size(),data,weights,wmean,wsd);
1674  }
1675
1676  /** \brief Compute the kurtosis of data with specified mean
1677  and standard deviation
1678
1679  This function produces the same results
1680  as <tt>gsl_stats_wkurtosis()</tt>.
1681  */
1682  template<class vec_t, class vec2_t>
1683  double wvector_kurtosis(size_t n, const vec_t &data,
1684  const vec2_t &weights) {
1685  double wmean=wvector_mean(n,data,weights);
1686  double wsd=wvector_stddev(n,data,weights,wmean);
1687  return wvector_kurtosis(n,data,weights,wmean,wsd);
1688  }
1689
1690  /** \brief Compute the kurtosis of data with specified mean
1691  and standard deviation
1692
1693  This function produces the same results
1694  as <tt>gsl_stats_wkurtosis()</tt>.
1695  */
1696  template<class vec_t, class vec2_t>
1697  double wvector_kurtosis(const vec_t &data, const vec2_t &weights) {
1698  return wvector_kurtosis<vec_t,vec2_t>(data,weights);
1699  }
1700  //@}
1701
1702 #ifndef DOXYGEN_NO_O2NS
1703 }
1704 #endif
1705
1706 #endif
double vector_lagk_autocorr(size_t n, const vec_t &data, size_t k, double mean)
Lag-k autocorrelation.
Definition: vec_stats.h:705
double vector_mean(size_t n, const vec_t &data)
Compute the mean of the first n elements of a vector.
Definition: vec_stats.h:55
double wvector_mean(size_t n, const vec_t &data, const vec2_t &weights)
Compute the mean of weighted data.
Definition: vec_stats.h:1186
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
double wvector_absdev(size_t n, const vec_t &data, const vec2_t &weights, double wmean)
Compute the absolute deviation of data about a specified mean.
Definition: vec_stats.h:1534
double wvector_variance_fmean(size_t n, const vec_t &data, const vec2_t &weights, double wmean)
Compute the variance of a weighted vector with a mean known in advance.
Definition: vec_stats.h:1269
double vector_correlation(size_t n, const vec_t &data1, const vec2_t &data2)
Pearson&#39;s correlation.
Definition: vec_stats.h:907
double vector_variance_fmean(size_t n, const vec_t &data, double mean)
Compute variance with specified mean known in advance.
Definition: vec_stats.h:87
invalid argument supplied by user
Definition: err_hnd.h:59
double vector_skew(size_t n, const vec_t &data, double mean, double stddev)
Skewness with specified mean and standard deviation.
Definition: vec_stats.h:428
double vector_kurtosis(size_t n, const vec_t &data, double mean, double stddev)
Kurtosis with specified mean and standard deviation.
Definition: vec_stats.h:518
double vector_absdev(size_t n, const vec_t &data, double mean)
Absolute deviation from the specified mean.
Definition: vec_stats.h:340
double wvector_kurtosis(size_t n, const vec_t &data, const vec2_t &weights, double wmean, double wsd)
Compute the kurtosis of data with specified mean and standard deviation.
Definition: vec_stats.h:1648
double vector_variance(size_t n, const vec_t &data, double mean)
Compute the variance with specified mean.
Definition: vec_stats.h:128
double vector_pvariance(size_t n1, const vec_t &data1, size_t n2, const vec2_t &data2)
The pooled variance of two vectors.
Definition: vec_stats.h:994
double vector_median_sorted(size_t n, const vec_t &data)
Return the median of sorted (ascending or descending) data.
Definition: vec_stats.h:1096
double vector_stddev(size_t n, const vec_t &data)
Standard deviation with specified mean.
Definition: vec_stats.h:253
double wvector_stddev_fmean(size_t n, const vec_t &data, const vec2_t &weights, double wmean)
Compute the standard deviation of a weighted vector with a mean known in advance. ...
Definition: vec_stats.h:1400
double wvector_covariance(size_t n, const vec_t &data1, const vec2_t &data2, const vec3_t &weights)
The weighted covariance of two vectors.
Definition: vec_stats.h:1363
double wvector_stddev(size_t n, const vec_t &data, const vec2_t &weights)
Compute the standard deviation of a weighted vector where mean is computed automatically.
Definition: vec_stats.h:1425
double wvector_variance(size_t n, const vec_t &data, const vec2_t &weights, double wmean)
Compute the variance of a weighted vector with specified mean.
Definition: vec_stats.h:1311
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
double wvector_sumsq(size_t n, const vec_t &data, const vec2_t &weights, double wmean)
Compute the weighted sum of squares of data about the specified weighted mean.
Definition: vec_stats.h:1477
double vector_covariance(size_t n, const vec_t &data1, const vec2_t &data2, double mean1, double mean2)
Compute the covariance of two vectors.
Definition: vec_stats.h:800
double vector_quantile_sorted(size_t n, const vec_t &data, const double f)
Quantile from sorted data (ascending only)
Definition: vec_stats.h:1044
double vector_lag1_autocorr(size_t n, const vec_t &data, double mean)
Lag-1 autocorrelation.
Definition: vec_stats.h:609
double vector_chi_squared(size_t n, const vec_t &obs, const vec2_t &exp, const vec3_t &err)
Compute the chi-squared statistic.
Definition: vec_stats.h:1141
double wvector_factor(size_t n, const vec_t &weights)
Compute a normalization factor for weighted data.
Definition: vec_stats.h:1228
double wvector_skew(size_t n, const vec_t &data, const vec2_t &weights, double wmean, double wsd)
Compute the skewness of data with specified mean and standard deviation.
Definition: vec_stats.h:1590
double vector_stddev_fmean(size_t n, const vec_t &data, double mean)
Standard deviation with specified mean known in advance.
Definition: vec_stats.h:215

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).