vec_stats.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, 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 mean, std. dev., and variance in src/other/vec_stats.h
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 
353  /// \name Vector absolute deviation, skewness, etc. in src/other/vec_stats.h
354  //@{
355  /** \brief Absolute deviation from the specified mean
356 
357  This function computes
358  \f[
359  \sum_i | x_i - \mu |
360  \f]
361  where the value of \f$ \mu \f$ is given in \c mean.
362 
363  This function produces the same results
364  as <tt>gsl_stats_absdev_m()</tt>.
365 
366  If \c n is zero, this function will return zero
367  without calling the error handler.
368  */
369  template<class vec_t> double vector_absdev(const vec_t &data,
370  double mean) {
371  return vector_absdev(data.size(),data,mean);
372  }
373 
374  /** \brief Absolute deviation from the computed mean
375 
376  This function computes
377  \f[
378  \sum_i | x_i - \mu |
379  \f]
380  where the value of \f$ \mu \f$ is mean as computed
381  from \ref vector_mean().
382 
383  This function produces the same results
384  as <tt>gsl_stats_absdev()</tt>.
385 
386  If \c n is zero, this function will return zero
387  without calling the error handler.
388  */
389  template<class vec_t>
390  double vector_absdev(size_t n, const vec_t &data) {
391  double mean=vector_mean<vec_t>(n,data);
392  return vector_absdev(n,data,mean);
393  }
394 
395  /** \brief Absolute deviation from the computed mean
396 
397  This function computes
398  \f[
399  \sum_i | x_i - \mu |
400  \f]
401  where the value of \f$ \mu \f$ is mean as computed
402  from \ref vector_mean().
403 
404  This function produces the same results
405  as <tt>gsl_stats_absdev()</tt>.
406 
407  If \c n is zero, this function will return zero
408  without calling the error handler.
409  */
410  template<class vec_t>
411  double vector_absdev(const vec_t &data) {
412  return vector_absdev(data.size(),data);
413  }
414 
415  /** \brief Skewness with specified mean and standard deviation
416 
417  This function computes
418  \f[
419  \frac{1}{N} \sum_i \left[
420  \frac{ \left(x_i - \mu \right)}{ \sigma }\right]^3
421  \f]
422  where the values of \f$ \mu \f$ and \f$ \sigma \f$
423  are given in \c mean and \c stddev.
424 
425  This function produces the same results
426  as <tt>gsl_stats_skew_m_sd()</tt>.
427 
428  If \c n is zero, this function will return zero
429  without calling the error handler.
430  */
431  template<class vec_t> double vector_skew(size_t n, const vec_t &data,
432  double mean, double stddev) {
433  long double skew=0.0;
434  for(size_t i=0;i<n;i++) {
435  long double x=(data[i]-mean)/stddev;
436  skew+=(x*x*x-skew)/(i+1);
437  }
438  return skew;
439  }
440 
441  /** \brief Skewness with specified mean and standard deviation
442 
443  This function computes
444  \f[
445  \frac{1}{N} \sum_i \left[
446  \frac{ \left(x_i - \mu \right)}{ \sigma }\right]^3
447  \f]
448  where the values of \f$ \mu \f$ and \f$ \sigma \f$
449  are given in \c mean and \c stddev.
450 
451  This function produces the same results
452  as <tt>gsl_stats_skew_m_sd()</tt>.
453 
454  If \c n is zero, this function will return zero
455  without calling the error handler.
456  */
457  template<class vec_t> double vector_skew(const vec_t &data,
458  double mean, double stddev) {
459  return vector_skew(data.size(),data,mean,stddev);
460  }
461 
462  /** \brief Skewness with computed mean and standard deviation
463 
464  This function computes
465  \f[
466  \frac{1}{N} \sum_i \left[
467  \frac{ \left(x_i - \mu \right)}{ \sigma }\right]^3
468  \f]
469  where the values of \f$ \mu \f$ and \f$ \sigma \f$
470  are computed using \ref vector_mean() and \ref vector_stddev().
471 
472  This function produces the same results
473  as <tt>gsl_stats_skew()</tt>.
474 
475  If \c n is zero, this function will return zero
476  without calling the error handler.
477  */
478  template<class vec_t> double vector_skew(size_t n, const vec_t &data) {
479  double mean=vector_mean<vec_t>(n,data);
480  double sd=vector_stddev<vec_t>(n,data,mean);
481  return vector_skew(n,data,mean,sd);
482  }
483 
484  /** \brief Skewness with computed mean and standard deviation
485 
486  This function computes
487  \f[
488  \frac{1}{N} \sum_i \left[
489  \frac{ \left(x_i - \mu \right)}{ \sigma }\right]^3
490  \f]
491  where the values of \f$ \mu \f$ and \f$ \sigma \f$
492  are computed using \ref vector_mean() and \ref vector_stddev().
493 
494  This function produces the same results
495  as <tt>gsl_stats_skew()</tt>.
496 
497  If \c n is zero, this function will return zero
498  without calling the error handler.
499  */
500  template<class vec_t> double vector_skew(const vec_t &data) {
501  return vector_skew(data.size(),data);
502  }
503 
504  /** \brief Kurtosis with specified mean and standard deviation
505 
506  This function computes
507  \f[
508  -3 + \frac{1}{N} \sum_i \left[
509  \frac{ \left(x_i - \mu \right)}{ \sigma }\right]^4
510  \f]
511  where the values of \f$ \mu \f$ and \f$ \sigma \f$
512  are given in \c mean and \c stddev.
513 
514  This function produces the same results
515  as <tt>gsl_stats_kurtosis_m_sd()</tt>.
516 
517  If \c n is zero, this function will return zero
518  without calling the error handler.
519  */
520  template<class vec_t>
521  double vector_kurtosis(size_t n, const vec_t &data, double mean,
522  double stddev) {
523  long double avg=0.0;
524  for(size_t i=0;i<n;i++) {
525  long double x=(data[i]-mean)/stddev;
526  avg+=(x*x*x*x-avg)/(i+1);
527  }
528  return avg-3.0;
529  }
530 
531  /** \brief Kurtosis with specified mean and standard deviation
532 
533  This function computes
534  \f[
535  -3 + \frac{1}{N} \sum_i \left[
536  \frac{ \left(x_i - \mu \right)}{ \sigma }\right]^4
537  \f]
538  where the values of \f$ \mu \f$ and \f$ \sigma \f$
539  are given in \c mean and \c stddev.
540 
541  This function produces the same results
542  as <tt>gsl_stats_kurtosis_m_sd()</tt>.
543 
544  If \c n is zero, this function will return zero
545  without calling the error handler.
546  */
547  template<class vec_t>
548  double vector_kurtosis(const vec_t &data, double mean,
549  double stddev) {
550  return vector_kurtosis(data.size(),data,mean,stddev);
551  }
552 
553  /** \brief Kurtosis with computed mean and standard deviation
554 
555  This function computes
556  \f[
557  -3 + \frac{1}{N} \sum_i \left[
558  \frac{ \left(x_i - \mu \right)}{ \sigma }\right]^4
559  \f]
560  where the values of \f$ \mu \f$ and \f$ \sigma \f$
561  are computed using \ref vector_mean() and \ref vector_stddev().
562 
563  This function produces the same results
564  as <tt>gsl_stats_kurtosis()</tt>.
565 
566  If \c n is zero, this function will return zero
567  without calling the error handler.
568  */
569  template<class vec_t> double vector_kurtosis(size_t n, const vec_t &data) {
570  double mean=vector_mean<vec_t>(n,data);
571  double sd=vector_stddev<vec_t>(n,data,mean);
572  return vector_kurtosis(n,data,mean,sd);
573  }
574 
575  /** \brief Kurtosis with computed mean and standard deviation
576 
577  This function computes
578  \f[
579  -3 + \frac{1}{N} \sum_i \left[
580  \frac{ \left(x_i - \mu \right)}{ \sigma }\right]^4
581  \f]
582  where the values of \f$ \mu \f$ and \f$ \sigma \f$
583  are computed using \ref vector_mean() and \ref vector_stddev().
584 
585  This function produces the same results
586  as <tt>gsl_stats_kurtosis()</tt>.
587 
588  If \c n is zero, this function will return zero
589  without calling the error handler.
590  */
591  template<class vec_t> double vector_kurtosis(const vec_t &data) {
592  return vector_kurtosis(data.size(),data);
593  }
594  //@}
595 
596  /// \name Other vector functions in src/other/vec_stats.h
597  //@{
598  /** \brief Compute the covariance of two vectors
599 
600  This function computes
601  \f[
602  \frac{1}{n-1} \sum_i \left(x_i - {\bar{x}}\right)
603  \left(y_i - {\bar{y}}\right)
604  \f]
605  where \f$ {\bar{x}} \f$ and \f$ {\bar{y}} \f$ are specified
606  in \c mean1 and \c mean2, respectively.
607 
608  This function produces the same results
609  as <tt>gsl_stats_covariance_m()</tt>.
610 
611  If \c n is zero, this function will return zero
612  without calling the error handler.
613  */
614  template<class vec_t, class vec2_t>
615  double vector_covariance(size_t n, const vec_t &data1, const vec2_t &data2,
616  double mean1, double mean2) {
617  double covar=0.0;
618  for(size_t i=0;i<n;i++) {
619  double delta1=(data1[i]-mean1);
620  double delta2=(data2[i]-mean2);
621  covar+=(delta1*delta2-covar)/(i+1);
622  }
623  return covar*n/(n-1);
624  }
625 
626  /** \brief Compute the covariance of two vectors
627 
628  This function computes
629  \f[
630  \frac{1}{n-1} \sum_i \left(x_i - {\bar{x}}\right)
631  \left(y_i - {\bar{y}}\right)
632  \f]
633  where \f$ {\bar{x}} \f$ and \f$ {\bar{y}} \f$ are specified
634  in \c mean1 and \c mean2, respectively.
635 
636  This function produces the same results
637  as <tt>gsl_stats_covariance_m()</tt>.
638 
639  If \c n is zero, this function will return zero
640  without calling the error handler.
641  */
642  template<class vec_t, class vec2_t>
643  double vector_covariance(const vec_t &data1, const vec2_t &data2,
644  double mean1, double mean2) {
645  return vector_covariance(data1.size(),data1,data2,mean1,mean2);
646  }
647 
648  /** \brief Compute the covariance of two vectors
649 
650  This function computes
651  \f[
652  \frac{1}{n-1} \sum_i \left(x_i - {\bar{x}}\right)
653  \left(y_i - {\bar{y}}\right)
654  \f]
655  where \f$ {\bar{x}} \f$ and \f$ {\bar{y}} \f$ are
656  the averages of \c data1 and \c data2 and are computed
657  automatically using \ref vector_mean().
658 
659  This function produces the same
660  results as <tt>gsl_stats_covariance()</tt>.
661 
662  If \c n is zero, this function will return zero
663  without calling the error handler.
664  */
665  template<class vec_t, class vec2_t>
666  double vector_covariance(size_t n, const vec_t &data1,
667  const vec2_t &data2) {
668  double covar=0.0;
669  double mean1=vector_mean<vec_t>(n,data1);
670  double mean2=vector_mean<vec_t>(n,data2);
671  for(size_t i=0;i<n;i++) {
672  long double delta1=(data1[i]-mean1);
673  long double delta2=(data2[i]-mean2);
674  covar+=(delta1*delta2-covar)/(i+1);
675  }
676  return covar*n/(n-1);
677  }
678 
679  /** \brief Compute the covariance of two vectors
680 
681  This function computes
682  \f[
683  \frac{1}{n-1} \sum_i \left(x_i - {\bar{x}}\right)
684  \left(y_i - {\bar{y}}\right)
685  \f]
686  where \f$ {\bar{x}} \f$ and \f$ {\bar{y}} \f$ are
687  the averages of \c data1 and \c data2 and are computed
688  automatically using \ref vector_mean().
689 
690  This function produces the same
691  results as <tt>gsl_stats_covariance()</tt>.
692 
693  If \c n is zero, this function will return zero
694  without calling the error handler.
695  */
696  template<class vec_t, class vec2_t>
697  double vector_covariance(const vec_t &data1,
698  const vec2_t &data2) {
699  return vector_covariance(data1.size(),data1,data2);
700  }
701 
702  /** \brief Pearson's correlation
703 
704  This function computes the Pearson correlation coefficient
705  between \c data1 and \c data2 .
706 
707  This function produces the same
708  results as <tt>gsl_stats_correlation()</tt>.
709 
710  \comment
711  r = cov(x, y) / (\Hat\sigma_x \Hat\sigma_y)
712  = {1/(n-1) \sum (x_i - \Hat x) (y_i - \Hat y)
713  \over
714  \sqrt{1/(n-1) \sum (x_i - \Hat x)^2} \sqrt{1/(n-1)
715  \sum (y_i - \Hat y)^2}
716  }
717  \endcomment
718 
719  If \c n is zero, this function will call the error handler.
720  */
721  template<class vec_t, class vec2_t>
722  double vector_correlation(size_t n, const vec_t &data1,
723  const vec2_t &data2) {
724  size_t i;
725 
726  if (n<1) {
727  O2SCL_ERR2("Cannot compute correlation with no elements",
728  " in vector_correlation().",exc_einval);
729  }
730 
731  double sum_xsq=0.0;
732  double sum_ysq=0.0;
733  double sum_cross=0.0;
734  double ratio;
735  double delta_x, delta_y;
736  double mean_x, mean_y;
737  double r;
738 
739  /*
740  * Compute:
741  * sum_xsq = Sum [ (x_i - mu_x)^2 ],
742  * sum_ysq = Sum [ (y_i - mu_y)^2 ] and
743  * sum_cross = Sum [ (x_i - mu_x) * (y_i - mu_y) ]
744  * using the above relation from Welford's paper
745  */
746 
747  mean_x=data1[0];
748  mean_y=data2[0];
749 
750  for (i=1; i < n; ++i) {
751  ratio=i / (i + 1.0);
752  delta_x=data1[i] - mean_x;
753  delta_y=data2[i] - mean_y;
754  sum_xsq += delta_x * delta_x * ratio;
755  sum_ysq += delta_y * delta_y * ratio;
756  sum_cross += delta_x * delta_y * ratio;
757  mean_x += delta_x / (i + 1.0);
758  mean_y += delta_y / (i + 1.0);
759  }
760 
761  r=sum_cross / (std::sqrt(sum_xsq) * std::sqrt(sum_ysq));
762 
763  return r;
764  }
765 
766  /** \brief Pearson's correlation
767 
768  This function computes the Pearson correlation coefficient
769  between \c data1 and \c data2 .
770 
771  This function produces the same
772  results as <tt>gsl_stats_correlation()</tt>.
773 
774  \comment
775  r = cov(x, y) / (\Hat\sigma_x \Hat\sigma_y)
776  = {1/(n-1) \sum (x_i - \Hat x) (y_i - \Hat y)
777  \over
778  \sqrt{1/(n-1) \sum (x_i - \Hat x)^2} \sqrt{1/(n-1)
779  \sum (y_i - \Hat y)^2}
780  }
781  \endcomment
782 
783  If \c n is zero, this function will call the error handler.
784  */
785  template<class vec_t, class vec2_t>
786  double vector_correlation(const vec_t &data1,
787  const vec2_t &data2) {
788  return vector_correlation(data1.size(),data1,data2);
789  }
790 
791  /** \brief The pooled variance of two vectors
792 
793  This function computes
794  \f[
795  s_{p}^2 = \frac{(n_1-1)s_1^2+(n_2-1)s_2^2}{n_1+n_2-2}
796  \f]
797  where \f$ n_i \f$ is the number of elements in vector \f$ i \f$
798  and \f$ s_i^2 \f$ is the variance of vector \f$ i \f$.
799 
800  From http://en.wikipedia.org/wiki/Pooled_variance, "Under the
801  assumption of equal population variances, the pooled sample
802  variance provides a higher precision estimate of variance than
803  the individual sample variances."
804 
805  This function produces the same
806  results as <tt>gsl_stats_pvariance()</tt>.
807  */
808  template<class vec_t, class vec2_t>
809  double vector_pvariance(size_t n1, const vec_t &data1,
810  size_t n2, const vec2_t &data2) {
811  double var1=vector_variance<vec_t>(n1,data1);
812  double var2=vector_variance<vec2_t>(n2,data2);
813  return (((n1-1)*var1)+((n2-1)*var2))/(n1+n2-2);
814  }
815 
816  /** \brief The pooled variance of two vectors
817 
818  This function computes
819  \f[
820  s_{p}^2 = \frac{(n_1-1)s_1^2+(n_2-1)s_2^2}{n_1+n_2-2}
821  \f]
822  where \f$ n_i \f$ is the number of elements in vector \f$ i \f$
823  and \f$ s_i^2 \f$ is the variance of vector \f$ i \f$.
824 
825  From http://en.wikipedia.org/wiki/Pooled_variance, "Under the
826  assumption of equal population variances, the pooled sample
827  variance provides a higher precision estimate of variance than
828  the individual sample variances."
829 
830  This function produces the same
831  results as <tt>gsl_stats_pvariance()</tt>.
832  */
833  template<class vec_t, class vec2_t>
834  double vector_pvariance(const vec_t &data1,
835  const vec2_t &data2) {
836  return vector_pvariance(data1.size(),data1,data2.size(),data2);
837  }
838 
839  /** \brief Quantile from sorted data (ascending only)
840 
841  This function returns the quantile \c f of data which
842  has already been sorted in ascending order. The quantile,
843  \f$ q \f$ , is
844  found by interpolation using
845  \f[
846  q = \left(1-\delta\right) x_i \delta x_{i+1}
847  \f]
848  where \f$ i = \mathrm{floor}[ (n-1)f ] \f$ and
849  \f$ \delta = (n-1)f -i \f$ .
850 
851  This function produces the same
852  results as <tt>gsl_stats_quantile_from_sorted_data()</tt>.
853 
854  No checks are made to ensure the data is sorted, or to ensure
855  that \f$ 0 \leq 0 \leq 1 \f$. If \c n is zero, this function
856  will return zero without calling the error handler.
857  */
858  template<class vec_t>
859  double vector_quantile_sorted(size_t n, const vec_t &data,
860  const double f) {
861 
862  double index=f*(n-1);
863  size_t lhs=((size_t)index);
864  double delta=index-lhs;
865  if (n==0) return 0.0;
866  if (lhs==n-1) return data[lhs];
867  return (1-delta)*data[lhs]+delta*data[lhs+1];
868  }
869 
870  /** \brief Quantile from sorted data (ascending only)
871 
872  This function returns the quantile \c f of data which
873  has already been sorted in ascending order. The quantile,
874  \f$ q \f$ , is
875  found by interpolation using
876  \f[
877  q = \left(1-\delta\right) x_i \delta x_{i+1}
878  \f]
879  where \f$ i = \mathrm{floor}[ (n-1)f ] \f$ and
880  \f$ \delta = (n-1)f -i \f$ .
881 
882  This function produces the same
883  results as <tt>gsl_stats_quantile_from_sorted_data()</tt>.
884 
885  No checks are made to ensure the data is sorted, or to ensure
886  that \f$ 0 \leq 0 \leq 1 \f$. If \c n is zero, this function
887  will return zero without calling the error handler.
888  */
889  template<class vec_t>
890  double vector_quantile_sorted(const vec_t &data, const double f) {
891  return vector_quantile_sorted<vec_t>(data.size(),data,f);
892  }
893 
894  /** \brief Return the median of sorted (ascending or descending) data
895 
896  This function returns the median of sorted data (either
897  ascending or descending), assuming the data has already been
898  sorted. When the data set has an odd number of elements, the
899  median is the value of the element at index \f$ (n-1)/2 \f$,
900  otherwise, the median is taken to be the average of the elements
901  at indices \f$ (n-1)/2 \f$ and \f$ n/2 \f$ .
902 
903  This function produces the same
904  results as <tt>gsl_stats_median_from_sorted_data()</tt>.
905 
906  No checks are made to ensure the data is sorted. If \c n is
907  zero, this function will return zero without calling the error
908  handler.
909  */
910  template<class vec_t>
911  double vector_median_sorted(size_t n, const vec_t &data) {
912 
913  if (n==0) return 0.0;
914 
915  size_t lhs=(n-1)/2;
916  size_t rhs=n/2;
917 
918  if (lhs==rhs) return data[lhs];
919 
920  return (data[lhs]+data[rhs])/2.0;
921  }
922 
923  /** \brief Return the median of sorted (ascending or descending) data
924 
925  This function returns the median of sorted data (either
926  ascending or descending), assuming the data has already been
927  sorted. When the data set has an odd number of elements, the
928  median is the value of the element at index \f$ (n-1)/2 \f$,
929  otherwise, the median is taken to be the average of the elements
930  at indices \f$ (n-1)/2 \f$ and \f$ n/2 \f$ .
931 
932  This function produces the same
933  results as <tt>gsl_stats_median_from_sorted_data()</tt>.
934 
935  No checks are made to ensure the data is sorted. If \c n is
936  zero, this function will return zero without calling the error
937  handler.
938  */
939  template<class vec_t>
940  double vector_median_sorted(const vec_t &data) {
941  return vector_median_sorted<vec_t>(data.size(),data);
942  }
943 
944  /** \brief Compute the chi-squared statistic
945 
946  This function computes
947  \f[
948  \sum_i \left( \frac{\mathrm{obs}_i - \mathrm{exp}_i}
949  {\mathrm{err}_i}\right)^2
950  \f]
951  where \f$ \mathrm{obs} \f$ are the observed values,
952  \f$ \mathrm{exp} \f$ are the expected values, and
953  \f$ \mathrm{err} \f$ are the errors.
954  */
955  template<class vec_t, class vec2_t, class vec3_t>
956  double vector_chi_squared(size_t n, const vec_t &obs, const vec2_t &exp,
957  const vec3_t &err) {
958  double chi2=0.0;
959  for(size_t i=0;i<n;i++) {
960  chi2+=pow((obs[i]-exp[i])/err[i],2.0);
961  }
962  return chi2;
963  }
964 
965  /** \brief Compute the chi-squared statistic
966 
967  This function computes
968  \f[
969  \sum_i \left( \frac{\mathrm{obs}_i - \mathrm{exp}_i}
970  {\mathrm{err}_i}\right)^2
971  \f]
972  where \f$ \mathrm{obs} \f$ are the observed values,
973  \f$ \mathrm{exp} \f$ are the expected values, and
974  \f$ \mathrm{err} \f$ are the errors.
975  */
976  template<class vec_t, class vec2_t, class vec3_t>
977  double vector_chi_squared(const vec_t &obs, const vec2_t &exp,
978  const vec3_t &err) {
979  return vector_chi_squared<vec_t,vec2_t,vec3_t>(obs.size(),obs,exp,err);
980  }
981 
982  /** \brief Optimal bin size using Scott's method for the
983  first <tt>n</tt> elements
984  */
985  template<class vec_t> double vector_bin_size_scott
986  (size_t n, const vec_t &v) {
987  if (n<=1) return 0.0;
988  double ret=3.5*vector_stddev(n,v)/cbrt(((double)n));
989  return ret;
990  }
991 
992  /** \brief Optimal bin size using Scott's method
993 
994  This function computes the optimal bin size \f$ \Delta_b \f$ of
995  a histogram using the expression
996  \f[
997  \Delta_b = \frac{3.5 \sigma}{n^{1/3}}
998  \f]
999 
1000  From \ref Scott79 .
1001 
1002  \note If <tt>n</tt> is less than or equal to 1, this
1003  function returns 0.0 without calling the error handler.
1004  */
1005  template<class vec_t> double vector_bin_size_scott
1006  (const vec_t &v) {
1007  return vector_bin_size_scott(v.size(),v);
1008  }
1009 
1010  /** \brief Obtain a quantile from a sorted vector
1011 
1012  This is a generic version of
1013  <tt>gsl_stats_quantile_from_sorted_data()</tt>.
1014 
1015  If <tt>f</tt> is less than 0 or greater than 1, the error
1016  handler is called. If <tt>n</tt> is zero, this function returns
1017  zero without calling the error handler.
1018  */
1019  template<class vec_t> double vector_sorted_quantile
1020  (size_t n, const vec_t &v, double f) {
1021 
1022  if (f<0.0 || f>1.0) {
1023  O2SCL_ERR("Invalid fraction for vector_sorted_quantile",
1025  }
1026 
1027  double index=f*(n-1);
1028  size_t lhs=(int)index;
1029  double delta=index-lhs;
1030  double result;
1031 
1032  if (n == 0) {
1033  return 0.0;
1034  }
1035 
1036  if (lhs == n-1) {
1037  result=v[lhs];
1038  } else {
1039  result=(1-delta)*v[lhs]+delta*v[(lhs+1)];
1040  }
1041 
1042  return result;
1043  }
1044 
1045  /** \brief Optimal bin size using the Freedman-Diaconis rule
1046  for the first <tt>n</tt> elements
1047  */
1048  template<class vec_t> double vector_bin_size_freedman
1049  (size_t n, vec_t &v) {
1050  vector_sort<vec_t,double>(n,v);
1051  double ret=2.0*(vector_sorted_quantile(n,v,0.75)-
1052  vector_sorted_quantile(n,v,0.25))/cbrt(((double)n));
1053  return ret;
1054  }
1055 
1056  /** \brief Optimal bin size using the Freedman-Diaconis rule
1057 
1058  This function computes the optimal bin size \f$ \Delta_b \f$ of
1059  a histogram using the expression
1060  \f[
1061  \Delta_b = \frac{2\left(q_{0.75}-q_{0.25}\right)}{n^{1/3}}
1062  \f]
1063  where \f$ q_{i} \f$ is the \f$ i \f$ quantile
1064  of the data (note this is quantile not quartile).
1065  This function sorts the vector in order to obtain
1066  the result.
1067 
1068  From \ref Freedman81 .
1069 
1070  \note If <tt>n</tt> is less than or equal to 1, this
1071  function returns 0.0 without calling the error handler.
1072  */
1073  template<class vec_t> double vector_bin_size_freedman
1074  (vec_t &v) {
1075  return vector_bin_size_freedman(v.size(),v);
1076  }
1077  //@}
1078 
1079  /// \name Weighted vector mean, std. dev., etc. in src/other/vec_stats.h
1080  //@{
1081  /** \brief Compute the mean of weighted data
1082 
1083  This function computes
1084  \f[
1085  \left( \sum_i w_i x_i \right) \left( \sum_i w_i \right)^{-1}
1086  \f]
1087 
1088  This function produces the same results
1089  as <tt>gsl_stats_wmean()</tt>.
1090 
1091  \comment
1092  M(n) = M(n-1) + (data[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
1093  W(n) = W(n-1) + w(n)
1094  \endcomment
1095  */
1096  template<class vec_t, class vec2_t>
1097  double wvector_mean(size_t n, const vec_t &data, const vec2_t &weights) {
1098 
1099  long double wmean=0.0;
1100  long double W=0.0;
1101  for(size_t i=0;i<n;i++) {
1102  double wi=weights[i];
1103  if (wi>0.0) {
1104  W+=wi;
1105  wmean+=(data[i]-wmean)*(wi/W);
1106  }
1107  }
1108 
1109  return wmean;
1110  }
1111 
1112  /** \brief Compute the mean of weighted data
1113 
1114  This function computes
1115  \f[
1116  \left( \sum_i w_i x_i \right) \left( \sum_i w_i \right)^{-1}
1117  \f]
1118 
1119  This function produces the same results
1120  as <tt>gsl_stats_wmean()</tt>.
1121 
1122  \comment
1123  M(n) = M(n-1) + (data[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
1124  W(n) = W(n-1) + w(n)
1125  \endcomment
1126  */
1127  template<class vec_t, class vec2_t>
1128  double wvector_mean(const vec_t &data, const vec2_t &weights) {
1129  return wvector_mean<vec_t,vec2_t>(data.size(),data,weights);
1130  }
1131 
1132  /** \brief Compute a normalization factor for weighted data
1133 
1134  This function is used internally in \ref wvector_variance(size_t
1135  n, vec_t &data, const vec2_t &weights, double wmean) and \ref
1136  wvector_stddev(size_t n, vec_t &data, const vec2_t &weights, double
1137  wmean) .
1138  */
1139  template<class vec_t> double wvector_factor(size_t n, const vec_t &weights) {
1140 
1141  long double a=0.0;
1142  long double b=0.0;
1143  long double factor;
1144  for(size_t i=0;i<n;i++) {
1145  double wi=weights[i];
1146  if (wi>0.0) {
1147  a+=wi;
1148  b+=wi*wi;
1149  }
1150  }
1151  factor=a*a/(a*a-b);
1152  return factor;
1153  }
1154 
1155  /** \brief Compute a normalization factor for weighted data
1156 
1157  This function is used internally in \ref wvector_variance(size_t
1158  n, vec_t &data, const vec2_t &weights, double wmean) and \ref
1159  wvector_stddev(size_t n, vec_t &data, const vec2_t &weights, double
1160  wmean) .
1161  */
1162  template<class vec_t> double wvector_factor(const vec_t &weights) {
1163  return wvector_factor<vec_t>(weights.size(),weights);
1164  }
1165 
1166  /** \brief Compute the variance of a weighted vector with a mean
1167  known in advance
1168 
1169  This function computes
1170  \f[
1171  \left[ \sum_i w_i \left(x_i-\mu\right)^2 \right]
1172  \left[ \sum_i w_i \right]^{-1}
1173  \f]
1174 
1175  This function produces the same results
1176  as <tt>gsl_stats_wvariance_with_fixed_mean()</tt>.
1177 
1178  */
1179  template<class vec_t, class vec2_t>
1180  double wvector_variance_fmean(size_t n, const vec_t &data,
1181  const vec2_t &weights, double wmean) {
1182  long double wvariance=0.0;
1183  long double W=0.0;
1184  for(size_t i=0;i<n;i++) {
1185  double wi=weights[i];
1186  if (wi>0.0) {
1187  const long double delta=data[i]-wmean;
1188  W+=wi;
1189  wvariance+=(delta*delta-wvariance)*(wi/W);
1190  }
1191  }
1192 
1193  return wvariance;
1194  }
1195 
1196  /** \brief Compute the variance of a weighted vector with a mean
1197  known in advance
1198 
1199  This function computes
1200  \f[
1201  \left[ \sum_i w_i \left(x_i-\mu\right)^2 \right]
1202  \left[ \sum_i w_i \right]^{-1}
1203  \f]
1204 
1205  This function produces the same results
1206  as <tt>gsl_stats_wvariance_with_fixed_mean()</tt>.
1207 
1208  */
1209  template<class vec_t, class vec2_t>
1210  double wvector_variance_fmean(const vec_t &data,
1211  const vec2_t &weights, double wmean) {
1212  return wvector_variance_fmean(data.size(),data,weights,wmean);
1213  }
1214 
1215  /** \brief Compute the variance of a weighted vector with
1216  specified mean
1217 
1218  This function produces the same results
1219  as <tt>gsl_stats_wvariance_m()</tt>.
1220  */
1221  template<class vec_t, class vec2_t>
1222  double wvector_variance(size_t n, const vec_t &data,
1223  const vec2_t &weights, double wmean) {
1224 
1225  const double variance=wvector_variance_fmean
1226  (n,data,weights,wmean);
1227  const double scale=wvector_factor(n,weights);
1228  const double wvar=scale*variance;
1229  return wvar;
1230  }
1231 
1232  /** \brief Compute the variance of a weighted vector with
1233  specified mean
1234 
1235  This function produces the same results
1236  as <tt>gsl_stats_wvariance_m()</tt>.
1237  */
1238  template<class vec_t, class vec2_t>
1239  double wvector_variance(const vec_t &data,
1240  const vec2_t &weights, double wmean) {
1241  return wvector_variance<vec_t,vec2_t>(data.size(),data,weights,wmean);
1242  }
1243 
1244  /** \brief Compute the variance of a weighted vector where mean
1245  is computed automatically
1246 
1247  This function produces the same results
1248  as <tt>gsl_stats_wvariance()</tt>.
1249  */
1250  template<class vec_t, class vec2_t>
1251  double wvector_variance(size_t n, const vec_t &data,
1252  const vec2_t &weights) {
1253 
1254  double wmean=wvector_mean(n,data,weights);
1255  return wvector_variance<vec_t,vec2_t>(n,data,weights,wmean);
1256  }
1257 
1258  /** \brief Compute the variance of a weighted vector where mean
1259  is computed automatically
1260 
1261  This function produces the same results
1262  as <tt>gsl_stats_wvariance()</tt>.
1263  */
1264  template<class vec_t, class vec2_t>
1265  double wvector_variance(const vec_t &data, const vec2_t &weights) {
1266  return wvector_variance(data.size(),data,weights);
1267  }
1268 
1269  /** \brief Compute the standard deviation of a weighted vector
1270  with a mean known in advance
1271 
1272  This function produces the same results
1273  as <tt>gsl_stats_wsd_with_fixed_mean()</tt>.
1274  */
1275  template<class vec_t, class vec2_t>
1276  double wvector_stddev_fmean(size_t n, const vec_t &data,
1277  const vec2_t &weights, double wmean) {
1278  return sqrt(wvector_variance_fmean(n,data,weights,wmean));
1279  }
1280 
1281  /** \brief Compute the standard deviation of a weighted vector
1282  with a mean known in advance
1283 
1284  This function produces the same results
1285  as <tt>gsl_stats_wsd_with_fixed_mean()</tt>.
1286  */
1287  template<class vec_t, class vec2_t>
1288  double wvector_stddev_fmean(const vec_t &data,
1289  const vec2_t &weights, double wmean) {
1290  return wvector_stddev_fmean<vec_t,vec2_t>
1291  (data.size(),data,weights,wmean);
1292  }
1293 
1294  /** \brief Compute the standard deviation of a weighted vector where mean
1295  is computed automatically
1296 
1297  This function produces the same results
1298  as <tt>gsl_stats_wsd()</tt>.
1299  */
1300  template<class vec_t, class vec2_t>
1301  double wvector_stddev(size_t n, const vec_t &data,
1302  const vec2_t &weights) {
1303  double wmean=wvector_mean(n,data,weights);
1304  return sqrt(wvector_variance(n,data,weights,wmean));
1305  }
1306 
1307  /** \brief Compute the standard deviation of a weighted vector where mean
1308  is computed automatically
1309 
1310  This function produces the same results
1311  as <tt>gsl_stats_wsd()</tt>.
1312  */
1313  template<class vec_t, class vec2_t>
1314  double wvector_stddev(const vec_t &data, const vec2_t &weights) {
1315  return wvector_stddev(data.size(),data,weights);
1316  }
1317 
1318  /** \brief Compute the standard deviation of a weighted vector with
1319  specified mean
1320 
1321  This function produces the same results
1322  as <tt>gsl_stats_wsd_m()</tt>.
1323  */
1324  template<class vec_t, class vec2_t>
1325  double wvector_stddev(size_t n, const vec_t &data,
1326  const vec2_t &weights, double wmean) {
1327  const double variance=wvector_variance_fmean
1328  (n,data,weights,wmean);
1329  const double scale=wvector_factor(n,weights);
1330  double wvar=scale*variance;
1331  return sqrt(wvar);
1332  }
1333 
1334  /** \brief Compute the standard deviation of a weighted vector with
1335  specified mean
1336 
1337  This function produces the same results
1338  as <tt>gsl_stats_wsd_m()</tt>.
1339  */
1340  template<class vec_t, class vec2_t>
1341  double wvector_stddev(const vec_t &data,
1342  const vec2_t &weights, double wmean) {
1343  return wvector_stddev<vec_t,vec2_t>(data.size(),data,weights,wmean);
1344  }
1345  //@}
1346 
1347  /// \name Other weighted vector functions in src/other/vec_stats.h
1348  //@{
1349  /** \brief The weighted covariance of two vectors
1350 
1351  \note Experimental
1352  */
1353  template<class vec_t, class vec2_t, class vec3_t>
1354  double wvector_covariance(size_t n, const vec_t &data1,
1355  const vec2_t &data2,
1356  const vec3_t &weights) {
1357  double mean1=wvector_mean(n,data1,weights);
1358  double mean2=wvector_mean(n,data2,weights);
1359  double covar=0.0;
1360  double W=0.0;
1361  for(size_t i=0;i<n;i++) {
1362  double wi=weights[i];
1363  if (wi>0.0) {
1364  W+=wi;
1365  double delta1=(data1[i]-mean1);
1366  double delta2=(data2[i]-mean2);
1367  covar+=(wi/W)*(delta1*delta2-covar);
1368  }
1369  }
1370  double scale=wvector_factor(n,weights);
1371  return covar*scale;
1372  }
1373 
1374  /** \brief The weighted covariance of two vectors
1375 
1376  \note Experimental
1377  */
1378  template<class vec_t, class vec2_t, class vec3_t>
1379  double wvector_covariance(const vec_t &data1, const vec2_t &data2,
1380  const vec3_t &weights) {
1381  return wvector_covariance<vec_t,vec2_t,vec3_t>
1382  (data1.size(),data1,data2,weights);
1383  }
1384 
1385  /** \brief Compute the weighted sum of squares of data about the
1386  specified weighted mean
1387 
1388  This function produces the same results
1389  as <tt>gsl_stats_wtss_m()</tt>.
1390  */
1391  template<class vec_t, class vec2_t>
1392  double wvector_sumsq(size_t n, const vec_t &data,
1393  const vec2_t &weights, double wmean) {
1394  long double wtss=0.0;
1395  for(size_t i=0;i<n;i++) {
1396  double wi=weights[i];
1397  if (wi>0.0) {
1398  const long double delta=data[i]-wmean;
1399  wtss+=wi*delta*delta;
1400  }
1401  }
1402 
1403  return wtss;
1404  }
1405 
1406  /** \brief Compute the weighted sum of squares of data about the
1407  specified weighted mean
1408 
1409  This function produces the same results
1410  as <tt>gsl_stats_wtss_m()</tt>.
1411  */
1412  template<class vec_t, class vec2_t>
1413  double wvector_sumsq(const vec_t &data,
1414  const vec2_t &weights, double wmean) {
1415  return wvector_sumsq<vec_t,vec2_t>(data.size(),data,weights,wmean);
1416  }
1417 
1418  /** \brief Compute the weighted sum of squares of data about the
1419  weighted mean
1420 
1421  This function produces the same results
1422  as <tt>gsl_stats_wtss()</tt>.
1423  */
1424  template<class vec_t, class vec2_t>
1425  double wvector_sumsq(size_t n, const vec_t &data,
1426  const vec2_t &weights) {
1427 
1428  double wmean=wvector_mean(n,data,weights);
1429  return wvector_sumsq(n,data,weights,wmean);
1430  }
1431 
1432  /** \brief Compute the weighted sum of squares of data about the
1433  weighted mean
1434 
1435  This function produces the same results
1436  as <tt>gsl_stats_wtss()</tt>.
1437  */
1438  template<class vec_t, class vec2_t>
1439  double wvector_sumsq(const vec_t &data, const vec2_t &weights) {
1440  return wvector_sumsq<vec_t,vec2_t>(data.size(),data,weights);
1441  }
1442 
1443  /** \brief Compute the absolute deviation of data about a specified mean
1444 
1445  This function produces the same results
1446  as <tt>gsl_stats_wabsdev_m()</tt>.
1447  */
1448  template<class vec_t, class vec2_t>
1449  double wvector_absdev(size_t n, const vec_t &data, const vec2_t &weights,
1450  double wmean) {
1451  long double wabsdev=0.0;
1452  long double W=0.0;
1453  for(size_t i=0;i<n;i++) {
1454  double wi=weights[i];
1455  if (wi>0.0) {
1456  const long double delta=fabs(data[i]-wmean);
1457  W+=wi;
1458  wabsdev+=(delta-wabsdev)*(wi/W);
1459  }
1460  }
1461  return wabsdev;
1462  }
1463 
1464  /** \brief Compute the absolute deviation of data about a specified mean
1465 
1466  This function produces the same results
1467  as <tt>gsl_stats_wabsdev_m()</tt>.
1468  */
1469  template<class vec_t, class vec2_t>
1470  double wvector_absdev(const vec_t &data, const vec2_t &weights,
1471  double wmean) {
1472  return wvector_absdev<vec_t,vec2_t>(data.size(),data,weights,wmean);
1473  }
1474 
1475  /** \brief Compute the absolute deviation of data about a specified mean
1476 
1477  This function produces the same results
1478  as <tt>gsl_stats_wabsdev()</tt>.
1479  */
1480  template<class vec_t, class vec2_t>
1481  double wvector_absdev(size_t n, const vec_t &data,
1482  const vec2_t &weights) {
1483 
1484  double wmean=wvector_mean(n,data,weights);
1485  return wvector_absdev(n,data,weights,wmean);
1486  }
1487 
1488  /** \brief Compute the absolute deviation of data about a specified mean
1489 
1490  This function produces the same results
1491  as <tt>gsl_stats_wabsdev()</tt>.
1492  */
1493  template<class vec_t, class vec2_t>
1494  double wvector_absdev(const vec_t &data, const vec2_t &weights) {
1495  return wvector_absdev<vec_t,vec2_t>(data.size(),data,weights);
1496  }
1497 
1498  /** \brief Compute the skewness of data with specified mean
1499  and standard deviation
1500 
1501  This function produces the same results
1502  as <tt>gsl_stats_wskew_m_sd()</tt>.
1503  */
1504  template<class vec_t, class vec2_t>
1505  double wvector_skew(size_t n, const vec_t &data, const vec2_t &weights,
1506  double wmean, double wsd) {
1507  long double wskew=0.0;
1508  long double W=0.0;
1509  for(size_t i=0;i<n;i++) {
1510  double wi=weights[i];
1511  if (wi>0.0) {
1512  const long double x=(data[i]-wmean)/wsd;
1513  W+=wi;
1514  wskew+=(x*x*x-wskew)*(wi/W);
1515  }
1516  }
1517  return wskew;
1518  }
1519 
1520  /** \brief Compute the skewness of data with specified mean
1521  and standard deviation
1522 
1523  This function produces the same results
1524  as <tt>gsl_stats_wskew_m_sd()</tt>.
1525  */
1526  template<class vec_t, class vec2_t>
1527  double wvector_skew(const vec_t &data, const vec2_t &weights,
1528  double wmean, double wsd) {
1529  return wvector_skew<vec_t,vec2_t>(data.size(),data,weights,wmean,wsd);
1530  }
1531 
1532  /** \brief Compute the skewness of data with specified mean
1533  and standard deviation
1534 
1535  This function produces the same results
1536  as <tt>gsl_stats_wskew()</tt>.
1537  */
1538  template<class vec_t, class vec2_t>
1539  double wvector_skew(size_t n, const vec_t &data, const vec2_t &weights) {
1540  double wmean=wvector_mean(n,data,weights);
1541  double wsd=wvector_stddev(n,data,weights,wmean);
1542  return wvector_skew(n,data,weights,wmean,wsd);
1543  }
1544 
1545  /** \brief Compute the skewness of data with specified mean
1546  and standard deviation
1547 
1548  This function produces the same results
1549  as <tt>gsl_stats_wskew()</tt>.
1550  */
1551  template<class vec_t, class vec2_t>
1552  double wvector_skew(const vec_t &data, const vec2_t &weights) {
1553  return wvector_skew<vec_t,vec2_t>(data.size(),data,weights);
1554  }
1555 
1556  /** \brief Compute the kurtosis of data with specified mean
1557  and standard deviation
1558 
1559  This function produces the same results
1560  as <tt>gsl_stats_wkurtosis_m_sd()</tt>.
1561  */
1562  template<class vec_t, class vec2_t>
1563  double wvector_kurtosis(size_t n, const vec_t &data, const vec2_t &weights,
1564  double wmean, double wsd) {
1565  long double wavg=0.0;
1566  long double W=0.0;
1567  for(size_t i=0;i<n;i++) {
1568  double wi=weights[i];
1569  if (wi>0.0) {
1570  const long double x=(data[i]-wmean)/wsd;
1571  W+=wi;
1572  wavg+=(x*x*x*x-wavg)*(wi/W);
1573  }
1574  }
1575  return wavg-3.0;
1576  }
1577 
1578  /** \brief Compute the kurtosis of data with specified mean
1579  and standard deviation
1580 
1581  This function produces the same results
1582  as <tt>gsl_stats_wkurtosis_m_sd()</tt>.
1583  */
1584  template<class vec_t, class vec2_t>
1585  double wvector_kurtosis(const vec_t &data, const vec2_t &weights,
1586  double wmean, double wsd) {
1587  return wvector_kurtosis<vec_t,vec2_t>
1588  (data.size(),data,weights,wmean,wsd);
1589  }
1590 
1591  /** \brief Compute the kurtosis of data with specified mean
1592  and standard deviation
1593 
1594  This function produces the same results
1595  as <tt>gsl_stats_wkurtosis()</tt>.
1596  */
1597  template<class vec_t, class vec2_t>
1598  double wvector_kurtosis(size_t n, const vec_t &data,
1599  const vec2_t &weights) {
1600  double wmean=wvector_mean(n,data,weights);
1601  double wsd=wvector_stddev(n,data,weights,wmean);
1602  return wvector_kurtosis(n,data,weights,wmean,wsd);
1603  }
1604 
1605  /** \brief Compute the kurtosis of data with specified mean
1606  and standard deviation
1607 
1608  This function produces the same results
1609  as <tt>gsl_stats_wkurtosis()</tt>.
1610  */
1611  template<class vec_t, class vec2_t>
1612  double wvector_kurtosis(const vec_t &data, const vec2_t &weights) {
1613  return wvector_kurtosis<vec_t,vec2_t>(data,weights);
1614  }
1615  //@}
1616 
1617  // This section has to appear after wvector_mean()
1618  /// \name Vector autocorrelation in src/other/vec_stats.h
1619  //@{
1620  /** \brief Lag-1 autocorrelation
1621 
1622  This function computes
1623  \f[
1624  \left[
1625  \sum_i \left(x_i - \mu\right) \left(x_{i-1} - \mu \right)
1626  \right] \left[
1627  \sum_i \left(x_i - \mu\right)^2
1628  \right]^{-1}
1629  \f]
1630 
1631  This function produces the same results
1632  as <tt>gsl_stats_lag1_autocorrelation_m()</tt>.
1633 
1634  If \c n is less than 2, this function will call the error handler.
1635  */
1636  template<class vec_t>
1637  double vector_lag1_autocorr(size_t n, const vec_t &data, double mean) {
1638 
1639  if (n<2) {
1640  O2SCL_ERR2("Cannot compute lag1 with less than 2 elements",
1641  " in vector_lag1_autocorr().",exc_einval);
1642  }
1643 
1644  long double q=0.0;
1645  long double v=(data[0]-mean)*(data[0]-mean);
1646  for(size_t i=1;i<n;i++) {
1647  long double delta0=data[i-1]-mean;
1648  long double delta1=data[i]-mean;
1649  q+=(delta0*delta1-q)/(i+1);
1650  v+=(delta1*delta1-v)/(i+1);
1651  }
1652 
1653  return q/v;
1654  }
1655 
1656  /** \brief Lag-1 autocorrelation
1657 
1658  This function computes
1659  \f[
1660  \left[
1661  \sum_i \left(x_i - \mu\right) \left(x_{i-1} - \mu \right)
1662  \right] \left[
1663  \sum_i \left(x_i - \mu\right)^2
1664  \right]^{-1}
1665  \f]
1666 
1667  This function produces the same results
1668  as <tt>gsl_stats_lag1_autocorrelation_m()</tt>.
1669 
1670  If \c n is less than 2, this function will call the error handler.
1671  */
1672  template<class vec_t>
1673  double vector_lag1_autocorr(const vec_t &data, double mean) {
1674  return vector_lag1_autocorr(data.size(),data,mean);
1675  }
1676 
1677  /** \brief Lag-1 autocorrelation
1678 
1679  This function computes
1680  \f[
1681  \left[
1682  \sum_i \left(x_i - \mu\right) \left(x_{i-1} - \mu \right)
1683  \right] \left[
1684  \sum_i \left(x_i - \mu\right)^2
1685  \right]^{-1}
1686  \f]
1687 
1688  This function produces the same results
1689  as <tt>gsl_stats_lag1_autocorrelation()</tt>.
1690 
1691  If \c n is less than 2, this function will call the error handler.
1692  */
1693  template<class vec_t> double vector_lag1_autocorr
1694  (size_t n, const vec_t &data) {
1695  double mean=vector_mean<vec_t>(n,data);
1696  return vector_lag1_autocorr(n,data,mean);
1697  }
1698 
1699  /** \brief Lag-1 autocorrelation
1700 
1701  This function computes
1702  \f[
1703  \left[
1704  \sum_i \left(x_i - \mu\right) \left(x_{i-1} - \mu \right)
1705  \right] \left[
1706  \sum_i \left(x_i - \mu\right)^2
1707  \right]^{-1}
1708  \f]
1709 
1710  This function produces the same results
1711  as <tt>gsl_stats_lag1_autocorrelation()</tt>.
1712 
1713  If \c n is less than 2, this function will call the error handler.
1714  */
1715  template<class vec_t> double vector_lag1_autocorr(const vec_t &data) {
1716  return vector_lag1_autocorr(data.size(),data);
1717  }
1718 
1719  /** \brief Lag-k autocorrelation
1720 
1721  This function computes
1722  \f[
1723  \left[
1724  \sum_i \left(x_i - \mu\right) \left(x_{i-k} - \mu \right)
1725  \right] \left[
1726  \sum_i \left(x_i - \mu\right)^2
1727  \right]^{-1}
1728  \f]
1729 
1730  If <tt>n<=k</tt>, this function will call the error handler.
1731  */
1732  template<class vec_t>
1733  double vector_lagk_autocorr(size_t n, const vec_t &data, size_t k,
1734  double mean) {
1735 
1736  if (n<=k) {
1737  O2SCL_ERR2("Not enough elements ",
1738  "in vector_lagk_autocorr().",exc_einval);
1739  }
1740 
1741  long double q=0.0, v=0.0;
1742  for(size_t i=0;i<k;i++) {
1743  v+=(data[i]-mean)*(data[i]-mean)/(i+1);
1744  }
1745  for(size_t i=k;i<n;i++) {
1746  long double delta0=data[i-k]-mean;
1747  long double delta1=data[i]-mean;
1748  q+=(delta0*delta1-q)/(i+1);
1749  v+=(delta1*delta1-v)/(i+1);
1750  }
1751  return q/v;
1752  }
1753 
1754  /** \brief Lag-k autocorrelation
1755 
1756  This function computes
1757  \f[
1758  \left[
1759  \sum_i \left(x_i - \mu\right) \left(x_{i-k} - \mu \right)
1760  \right] \left[
1761  \sum_i \left(x_i - \mu\right)^2
1762  \right]^{-1}
1763  \f]
1764 
1765  If <tt>n<=k</tt>, this function will call the error handler.
1766  */
1767  template<class vec_t>
1768  double vector_lagk_autocorr(const vec_t &data, size_t k,
1769  double mean) {
1770  return vector_lagk_autocorr(data.size(),k,mean);
1771  }
1772 
1773  /** \brief Lag-k autocorrelation
1774 
1775  This function computes
1776  \f[
1777  \left[
1778  \sum_i \left(x_i - \mu\right) \left(x_{i-k} - \mu \right)
1779  \right] \left[
1780  \sum_i \left(x_i - \mu\right)^2
1781  \right]^{-1}
1782  \f]
1783 
1784  If <tt>n<=k</tt>, this function will call the error handler.
1785  */
1786  template<class vec_t> double vector_lagk_autocorr
1787  (size_t n, const vec_t &data, size_t k) {
1788  double mean=vector_mean<vec_t>(n,data);
1789  return vector_lagk_autocorr(n,data,k,mean);
1790  }
1791 
1792  /** \brief Lag-k autocorrelation
1793 
1794  This function computes
1795  \f[
1796  \left[
1797  \sum_i \left(x_i - \mu\right) \left(x_{i-k} - \mu \right)
1798  \right] \left[
1799  \sum_i \left(x_i - \mu\right)^2
1800  \right]^{-1}
1801  \f]
1802 
1803  If <tt>n<=k</tt>, this function will call the error handler.
1804  */
1805  template<class vec_t> double vector_lagk_autocorr
1806  (const vec_t &data, size_t k) {
1807  return vector_lagk_autocorr(data.size(),data,k);
1808  }
1809 
1810  /** \brief Construct an autocorrelation vector
1811 
1812  This constructs a vector \c ac_vec for which the kth entry
1813  stores the lag-k autocorrelation. This function chooses \f$
1814  k_{\mathrm{max}} =n/2 \f$ where \f$ n \f$ is the length of the
1815  \c data vector. The vector \c ac_vec is resized to accomodate
1816  exactly \f$ k_{\mathrm{max}} \f$ values, from 0 to
1817  \f$ k_{\mathrm{max}}-1 \f$.
1818  */
1819  template<class vec_t, class resize_vec_t> void vector_autocorr_vector
1820  (const vec_t &data, resize_vec_t &ac_vec) {
1821  size_t kmax=data.size()/2;
1822  double mean=vector_mean(data);
1823  ac_vec.resize(kmax);
1824  ac_vec[0]=1.0;
1825  for(size_t k=1;k<kmax;k++) {
1826  ac_vec[k]=vector_lagk_autocorr(data.size(),data,k,mean);
1827  }
1828  return;
1829  }
1830 
1831  /** \brief Use the Goodman method to compute the
1832  autocorrelation length
1833 
1834  Representing the lag-k correlation coefficient by
1835  \f$ \hat{C}(k) \f$, Goodman defines
1836  \f[
1837  \hat{\tau}(M) = 1 + 2 \sum_{s=1}^{M} \frac{\hat{C}(k)}{\hat{C}(0)}
1838  \f]
1839  Then the autocorrelation length is the value of
1840  \f$ \hat{\tau}(M) \f$ for which
1841  \f[
1842  5 \hat{\tau}(M)/M \leq 1
1843  \f]
1844 
1845  This function computes the value of \f$ 5 \hat{\tau}(M)/M \f$
1846  and stores it in the <tt>five_tau_over_M</tt> vector and then
1847  returns the first value of \f$ M \f$ for which the vector is
1848  less than or equal to 1.0. If this function returns 0, then \f$
1849  5 \hat{\tau}(M)/M \f$ is greater than 1.0 for all \f$ M \f$, and
1850  this can be a sign that the autocorrelation length is too long
1851  to accurately resolve.
1852 
1853  On completion, the vector \c five_tau_over_m will have
1854  one less element than the vector \c ac_vec .
1855  */
1856  template<class vec_t, class resize_vec_t> size_t vector_autocorr_tau
1857  (const vec_t &ac_vec, resize_vec_t &five_tau_over_M) {
1858  five_tau_over_M.resize(0);
1859  size_t len=0;
1860  bool len_set=false;
1861  for (size_t M=1;M<ac_vec.size();M++) {
1862  double sum=0.0;
1863  for(size_t s=1;s<=M;s++) {
1864  sum+=ac_vec[s];
1865  }
1866  double val=(1.0+2.0*sum)/((double)M)*5.0;
1867  if (len_set==false && val<=1.0) {
1868  len=M;
1869  len_set=true;
1870  }
1871  five_tau_over_M.push_back(val);
1872  }
1873  return len;
1874  }
1875 
1876  /** \brief Lag-k autocorrelation for the first
1877  \c n elements with a vector multiplier given
1878  the mean
1879  */
1880  template<class vec_t, class vec2_t>
1881  double vector_lagk_autocorr_mult(size_t n, const vec_t &data,
1882  const vec2_t &mult, size_t k,
1883  double mean) {
1884 
1885  size_t n2=0;
1886  for(size_t i=0;i<n;i++) {
1887  size_t m=((size_t)(mult[i]*(1.0+1.0e-10)));
1888  if (m==0) {
1889  O2SCL_ERR2("Mult vector is zero ",
1890  "in vector_lagk_autocorr_mult().",exc_einval);
1891  }
1892  n2+=m;
1893  }
1894 
1895  if (n2<=k) {
1896  O2SCL_ERR2("Not enough elements ",
1897  "in vector_lagk_autocorr_mult().",exc_einval);
1898  }
1899 
1900  long double q=0.0, v=0.0;
1901  size_t im=0, ix=0, im2=0, ix2=0;
1902  for(size_t i=0;i<k;i++) {
1903  v+=(data[ix]-mean)*(data[ix]-mean)/(i+1);
1904  im++;
1905  if (im>=((size_t)(mult[ix]*(1.0+1.0e-10)))) {
1906  im=0;
1907  ix++;
1908  }
1909  }
1910  for(size_t i=k;i<n2;i++) {
1911  long double delta0=data[ix2]-mean;
1912  long double delta1=data[ix]-mean;
1913  q+=(delta0*delta1-q)/(i+1);
1914  v+=(delta1*delta1-v)/(i+1);
1915  im++;
1916  if (im>=((size_t)(mult[ix]*(1.0+1.0e-10)))) {
1917  im=0;
1918  ix++;
1919  }
1920  im2++;
1921  if (im2>=((size_t)(mult[ix2]*(1.0+1.0e-10)))) {
1922  im2=0;
1923  ix2++;
1924  }
1925  }
1926  return q/v;
1927  }
1928 
1929  /** \brief Lag-k autocorrelation for the first
1930  \c n elements with a vector multiplier
1931  */
1932  template<class vec_t, class vec2_t>
1933  double vector_lagk_autocorr_mult(size_t n, const vec_t &data,
1934  const vec2_t &mult, size_t k) {
1935  double mean=wvector_mean(n,mult,data);
1936  return vector_lagk_autocorr_mult(n,data,mult,k,mean);
1937  }
1938 
1939  /** \brief Lag-k autocorrelation with a vector multiplier
1940  given the mean
1941  */
1942  template<class vec_t, class vec2_t>
1943  double vector_lagk_autocorr_mult(const vec_t &data,
1944  const vec2_t &mult, size_t k,
1945  double mean) {
1946  return vector_lagk_autocorr_mult(data.size(),data,mult,k,mean);
1947  }
1948 
1949  /** \brief Lag-k autocorrelation with a vector multiplier
1950  */
1951  template<class vec_t, class vec2_t>
1952  double vector_lagk_autocorr_mult(const vec_t &data,
1953  const vec2_t &mult, size_t k) {
1954  return vector_lagk_autocorr_mult(data.size(),data,mult,k);
1955  }
1956 
1957  /** \brief Construct an autocorrelation vector using a multiplier
1958  using the first \c n2 elements of vectors \c data and \c mult
1959  */
1960  template<class vec_t, class vec2_t, class resize_vec_t>
1962  (size_t n2, const vec_t &data, const vec2_t &mult, resize_vec_t &ac_vec) {
1963 
1964  size_t n=0;
1965  for(size_t i=0;i<n2;i++) {
1966  n+=((size_t)(mult[i]*(1.0+1.0e-10)));
1967  }
1968 
1969  size_t kmax=n/2;
1970  double mean=wvector_mean(data,mult);
1971  ac_vec.resize(kmax);
1972  ac_vec[0]=1.0;
1973  for(size_t k=1;k<kmax;k++) {
1974  ac_vec[k]=vector_lagk_autocorr_mult(n2,data,mult,k,mean);
1975  }
1976  return;
1977  }
1978 
1979  /** \brief Construct an autocorrelation vector using a multiplier
1980  */
1981  template<class vec_t, class vec2_t, class resize_vec_t>
1983  (const vec_t &data, const vec2_t &mult, resize_vec_t &ac_vec) {
1984 
1985  vector_autocorr_vector_mult(data.size(),data,mult,ac_vec);
1986 
1987  return;
1988  }
1989  //@}
1990 
1991  /// \name Convert a vector to bin edges in src/other/vec_stats.h
1992  //@{
1993  /** \brief Take a vector of data and convert it to a vector
1994  of bin edges automatically adjusting for increasing or
1995  decreasing and linear or logarithmic spacing
1996  */
1997  template<class vec_t, class vec2_t>
1998  void vector_to_bins(const vec_t &v_grid, vec2_t &v_bins,
1999  int verbose=1) {
2000 
2001  size_t n=v_grid.size();
2002  v_bins.resize(n+1);
2003 
2004  // Vector of differences
2005  std::vector<double> diffs;
2006 
2007  // Compute quality factor for linear bins (smaller number is
2008  // better)
2009  vector_diffs(v_grid,diffs);
2010  double mean=vector_mean(diffs);
2011  double std=vector_stddev(diffs);
2012  double qual_lin=std/mean;
2013 
2014  // Compute quality factor for logarithmic bins (smaller number is
2015  // better)
2016  std::vector<double> logs;
2017  for(size_t i=0;i<n;i++) logs[i]=log(v_grid[i]);
2018  vector_diffs(logs,diffs);
2019  mean=vector_mean(diffs);
2020  std=vector_stddev(diffs);
2021  double qual_log=std/mean;
2022 
2023  if (qual_log<qual_lin) {
2024  if (verbose>0) {
2025  std::cout << "Auto-detected log mode in o2scl::vector_to_bins()."
2026  << std::endl;
2027  }
2028  if (logs[1]>logs[0]) {
2029  // Increasing, log
2030  v_bins[0]=exp(logs[0]-(logs[1]-logs[0])/2.0);
2031  v_bins[n]=exp(logs[n-1]+(logs[n-1]-logs[n-2])/2.0);
2032  } else {
2033  // Decreasing, log
2034  v_bins[0]=exp(logs[0]+(logs[0]-logs[1])/2.0);
2035  v_bins[n]=exp(logs[n-1]-(logs[n-2]-logs[n-1])/2.0);
2036  }
2037  for(size_t i=1;i<n-1;i++) {
2038  v_bins[i]=exp((logs[i-1]+logs[i])/2.0);
2039  }
2040  } else {
2041  if (v_grid[1]>v_grid[0]) {
2042  // Increasing, linear
2043  v_bins[0]=v_grid[0]-(v_grid[1]-v_grid[0])/2.0;
2044  v_bins[n]=v_grid[n-1]+(v_grid[n-1]-v_grid[n-2])/2.0;
2045  } else {
2046  // Decreasing, linear
2047  v_bins[0]=v_grid[0]+(v_grid[0]-v_grid[1])/2.0;
2048  v_bins[n]=v_grid[n-1]-(v_grid[n-2]-v_grid[n-1])/2.0;
2049  }
2050  for(size_t i=1;i<n-1;i++) {
2051  v_bins[i]=(v_grid[i-1]+v_grid[i])/2.0;
2052  }
2053  }
2054 
2055  return;
2056  }
2057  //@}
2058 
2059 #ifndef DOXYGEN_NO_O2NS
2060 }
2061 #endif
2062 
2063 #endif
o2scl::wvector_mean
double wvector_mean(size_t n, const vec_t &data, const vec2_t &weights)
Compute the mean of weighted data.
Definition: vec_stats.h:1097
o2scl::vector_bin_size_scott
double vector_bin_size_scott(size_t n, const vec_t &v)
Optimal bin size using Scott's method for the first n elements.
Definition: vec_stats.h:986
o2scl::vector_variance_fmean
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
o2scl::wvector_variance_fmean
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:1180
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
o2scl::vector_to_bins
void vector_to_bins(const vec_t &v_grid, vec2_t &v_bins, int verbose=1)
Take a vector of data and convert it to a vector of bin edges automatically adjusting for increasing ...
Definition: vec_stats.h:1998
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl::wvector_covariance
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:1354
o2scl::vector_mean
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
o2scl::vector_skew
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:431
o2scl::wvector_stddev_fmean
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:1276
o2scl::vector_lagk_autocorr_mult
double vector_lagk_autocorr_mult(size_t n, const vec_t &data, const vec2_t &mult, size_t k, double mean)
Lag-k autocorrelation for the first n elements with a vector multiplier given the mean.
Definition: vec_stats.h:1881
o2scl::vector_lag1_autocorr
double vector_lag1_autocorr(size_t n, const vec_t &data, double mean)
Lag-1 autocorrelation.
Definition: vec_stats.h:1637
o2scl::vector_bin_size_freedman
double vector_bin_size_freedman(size_t n, vec_t &v)
Optimal bin size using the Freedman-Diaconis rule for the first n elements
Definition: vec_stats.h:1049
o2scl::vector_autocorr_vector
void vector_autocorr_vector(const vec_t &data, resize_vec_t &ac_vec)
Construct an autocorrelation vector.
Definition: vec_stats.h:1820
o2scl::vector_variance
double vector_variance(size_t n, const vec_t &data, double mean)
Compute the variance with specified mean.
Definition: vec_stats.h:128
o2scl::vector_stddev_fmean
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
o2scl::wvector_sumsq
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:1392
o2scl::wvector_skew
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:1505
o2scl::vector_stddev
double vector_stddev(size_t n, const vec_t &data)
Standard deviation with specified mean.
Definition: vec_stats.h:253
o2scl::vector_median_sorted
double vector_median_sorted(size_t n, const vec_t &data)
Return the median of sorted (ascending or descending) data.
Definition: vec_stats.h:911
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::wvector_stddev
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:1301
o2scl::vector_chi_squared
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:956
o2scl::vector_correlation
double vector_correlation(size_t n, const vec_t &data1, const vec2_t &data2)
Pearson's correlation.
Definition: vec_stats.h:722
o2scl::vector_sorted_quantile
double vector_sorted_quantile(size_t n, const vec_t &v, double f)
Obtain a quantile from a sorted vector.
Definition: vec_stats.h:1020
o2scl::wvector_absdev
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:1449
o2scl::vector_covariance
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:615
o2scl::wvector_kurtosis
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:1563
o2scl::vector_autocorr_vector_mult
void vector_autocorr_vector_mult(size_t n2, const vec_t &data, const vec2_t &mult, resize_vec_t &ac_vec)
Construct an autocorrelation vector using a multiplier using the first n2 elements of vectors data an...
Definition: vec_stats.h:1962
o2scl::vector_absdev
double vector_absdev(size_t n, const vec_t &data, double mean)
Absolute deviation from the specified mean.
Definition: vec_stats.h:340
o2scl::vector_kurtosis
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:521
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::vector_autocorr_tau
size_t vector_autocorr_tau(const vec_t &ac_vec, resize_vec_t &five_tau_over_M)
Use the Goodman method to compute the autocorrelation length.
Definition: vec_stats.h:1857
o2scl::vector_lagk_autocorr
double vector_lagk_autocorr(size_t n, const vec_t &data, size_t k, double mean)
Lag-k autocorrelation.
Definition: vec_stats.h:1733
o2scl::wvector_factor
double wvector_factor(size_t n, const vec_t &weights)
Compute a normalization factor for weighted data.
Definition: vec_stats.h:1139
o2scl::vector_diffs
void vector_diffs(const vec_t &v_data, rvec_t &v_diffs)
Create a new vector containing the differences between adjacent entries.
Definition: vector.h:2193
o2scl::vector_pvariance
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:809
o2scl::vector_quantile_sorted
double vector_quantile_sorted(size_t n, const vec_t &data, const double f)
Quantile from sorted data (ascending only)
Definition: vec_stats.h:859
o2scl::wvector_variance
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:1222

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