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
9  it under the terms of the GNU General Public License as published by
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
1256  known in advance
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
1286  known in advance
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).