vec_stats.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2018, 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  v+=(data[i]-mean)*(data[i]-mean)/(i+1);
716  }
717  for(size_t i=k;i<n;i++) {
718  long double delta0=data[i-k]-mean;
719  long double delta1=data[i]-mean;
720  q+=(delta0*delta1-q)/(i+1);
721  v+=(delta1*delta1-v)/(i+1);
722  }
723  return q/v;
724  }
725 
726  /** \brief Lag-k autocorrelation
727 
728  This function computes
729  \f[
730  \left[
731  \sum_i \left(x_i - \mu\right) \left(x_{i-k} - \mu \right)
732  \right] \left[
733  \sum_i \left(x_i - \mu\right)^2
734  \right]^{-1}
735  \f]
736 
737  If <tt>n<=k</tt>, this function will call the error handler.
738  */
739  template<class vec_t>
740  double vector_lagk_autocorr(const vec_t &data, size_t k,
741  double mean) {
742  return vector_lagk_autocorr(data.size(),k,mean);
743  }
744 
745  /** \brief Lag-k autocorrelation
746 
747  This function computes
748  \f[
749  \left[
750  \sum_i \left(x_i - \mu\right) \left(x_{i-k} - \mu \right)
751  \right] \left[
752  \sum_i \left(x_i - \mu\right)^2
753  \right]^{-1}
754  \f]
755 
756  If <tt>n<=k</tt>, this function will call the error handler.
757  */
758  template<class vec_t> double vector_lagk_autocorr
759  (size_t n, const vec_t &data, size_t k) {
760  double mean=vector_mean<vec_t>(n,data);
761  return vector_lagk_autocorr(n,data,k,mean);
762  }
763 
764  /** \brief Lag-k autocorrelation
765 
766  This function computes
767  \f[
768  \left[
769  \sum_i \left(x_i - \mu\right) \left(x_{i-k} - \mu \right)
770  \right] \left[
771  \sum_i \left(x_i - \mu\right)^2
772  \right]^{-1}
773  \f]
774 
775  If <tt>n<=k</tt>, this function will call the error handler.
776  */
777  template<class vec_t> double vector_lagk_autocorr
778  (const vec_t &data, size_t k) {
779  return vector_lagk_autocorr(data.size(),data,k);
780  }
781 
782  /** \brief Construct an autocorrelation vector
783 
784  This constructs a vector \c ac_vec for which
785  the kth entry stores the lag-k autocorrelation.
786  */
787  template<class vec_t, class resize_vec_t> void vector_autocorr_vector
788  (const vec_t &data, resize_vec_t &ac_vec) {
789  size_t kmax=data.size()/2;
790  double mean=vector_mean(data);
791  ac_vec.resize(kmax);
792  ac_vec[0]=1.0;
793  for(size_t k=1;k<kmax;k++) {
794  ac_vec[k]=vector_lagk_autocorr(data.size(),data,k,mean);
795  }
796  return;
797  }
798 
799  /** \brief Use the Goodman method to compute the
800  autocorrelation length
801 
802  Representing the lag-k correlation coefficient by
803  \f$ \hat{C}(k) \f$, Goodman defines
804  \f[
805  \hat{\tau}(M) = 1 + 2 \sum_{s=1}^{M} \frac{\hat{C}(k)}{\hat{C}(0)}
806  \f]
807  Then the autocorrelation length is the value of
808  \f$ \hat{\tau}(M) \f$ for which
809  \f[
810  5 \hat{\tau}(M)/M \leq 1
811  \f]
812 
813  This function computes the value of \f$ 5 \hat{\tau}(M)/M \f$
814  and stores it in the <tt>five_tau_over_M</tt> vector and then
815  returns the first value of \f$ M \f$ for which the vector is
816  less than or equal to 1.0. If this function returns 0, then \f$
817  5 \hat{\tau}(M)/M \f$ is greater than 1.0 for all \f$ M \f$, and
818  this can be a sign that the autocorrelation length is too long
819  to accurately resolve.
820 
821  On completion, the vector \c five_tau_over_m will have
822  one less element than the vector \c ac_vec .
823  */
824  template<class vec_t, class resize_vec_t> size_t vector_autocorr_tau
825  (const vec_t &ac_vec, resize_vec_t &five_tau_over_M) {
826  five_tau_over_M.resize(0);
827  size_t len=0;
828  bool len_set=false;
829  for (size_t M=1;M<ac_vec.size();M++) {
830  double sum=0.0;
831  for(size_t s=1;s<=M;s++) {
832  sum+=ac_vec[s];
833  }
834  double val=(1.0+2.0*sum)/((double)M)*5.0;
835  if (len_set==false && val<=1.0) {
836  len=M;
837  len_set=true;
838  }
839  five_tau_over_M.push_back(val);
840  }
841  return len;
842  }
843 
844  /** \brief Compute the covariance of two vectors
845 
846  This function computes
847  \f[
848  \frac{1}{n-1} \sum_i \left(x_i - {\bar{x}}\right)
849  \left(y_i - {\bar{y}}\right)
850  \f]
851  where \f$ {\bar{x}} \f$ and \f$ {\bar{y}} \f$ are specified
852  in \c mean1 and \c mean2, respectively.
853 
854  This function produces the same results
855  as <tt>gsl_stats_covariance_m()</tt>.
856 
857  If \c n is zero, this function will return zero
858  without calling the error handler.
859  */
860  template<class vec_t, class vec2_t>
861  double vector_covariance(size_t n, const vec_t &data1, const vec2_t &data2,
862  double mean1, double mean2) {
863  double covar=0.0;
864  for(size_t i=0;i<n;i++) {
865  double delta1=(data1[i]-mean1);
866  double delta2=(data2[i]-mean2);
867  covar+=(delta1*delta2-covar)/(i+1);
868  }
869  return covar*n/(n-1);
870  }
871 
872  /** \brief Compute the covariance of two vectors
873 
874  This function computes
875  \f[
876  \frac{1}{n-1} \sum_i \left(x_i - {\bar{x}}\right)
877  \left(y_i - {\bar{y}}\right)
878  \f]
879  where \f$ {\bar{x}} \f$ and \f$ {\bar{y}} \f$ are specified
880  in \c mean1 and \c mean2, respectively.
881 
882  This function produces the same results
883  as <tt>gsl_stats_covariance_m()</tt>.
884 
885  If \c n is zero, this function will return zero
886  without calling the error handler.
887  */
888  template<class vec_t, class vec2_t>
889  double vector_covariance(const vec_t &data1, const vec2_t &data2,
890  double mean1, double mean2) {
891  return vector_covariance(data1.size(),data1,data2,mean1,mean2);
892  }
893 
894  /** \brief Compute the covariance of two vectors
895 
896  This function computes
897  \f[
898  \frac{1}{n-1} \sum_i \left(x_i - {\bar{x}}\right)
899  \left(y_i - {\bar{y}}\right)
900  \f]
901  where \f$ {\bar{x}} \f$ and \f$ {\bar{y}} \f$ are
902  the averages of \c data1 and \c data2 and are computed
903  automatically using \ref vector_mean().
904 
905  This function produces the same
906  results as <tt>gsl_stats_covariance()</tt>.
907 
908  If \c n is zero, this function will return zero
909  without calling the error handler.
910  */
911  template<class vec_t, class vec2_t>
912  double vector_covariance(size_t n, const vec_t &data1,
913  const vec2_t &data2) {
914  double covar=0.0;
915  double mean1=vector_mean<vec_t>(n,data1);
916  double mean2=vector_mean<vec_t>(n,data2);
917  for(size_t i=0;i<n;i++) {
918  long double delta1=(data1[i]-mean1);
919  long double delta2=(data2[i]-mean2);
920  covar+=(delta1*delta2-covar)/(i+1);
921  }
922  return covar*n/(n-1);
923  }
924 
925  /** \brief Compute the covariance of two vectors
926 
927  This function computes
928  \f[
929  \frac{1}{n-1} \sum_i \left(x_i - {\bar{x}}\right)
930  \left(y_i - {\bar{y}}\right)
931  \f]
932  where \f$ {\bar{x}} \f$ and \f$ {\bar{y}} \f$ are
933  the averages of \c data1 and \c data2 and are computed
934  automatically using \ref vector_mean().
935 
936  This function produces the same
937  results as <tt>gsl_stats_covariance()</tt>.
938 
939  If \c n is zero, this function will return zero
940  without calling the error handler.
941  */
942  template<class vec_t, class vec2_t>
943  double vector_covariance(const vec_t &data1,
944  const vec2_t &data2) {
945  return vector_covariance(data1.size(),data1,data2);
946  }
947 
948  /** \brief Pearson's correlation
949 
950  This function computes the Pearson correlation coefficient
951  between \c data1 and \c data2 .
952 
953  This function produces the same
954  results as <tt>gsl_stats_correlation()</tt>.
955 
956  \comment
957  r = cov(x, y) / (\Hat\sigma_x \Hat\sigma_y)
958  = {1/(n-1) \sum (x_i - \Hat x) (y_i - \Hat y)
959  \over
960  \sqrt{1/(n-1) \sum (x_i - \Hat x)^2} \sqrt{1/(n-1)
961  \sum (y_i - \Hat y)^2}
962  }
963  \endcomment
964 
965  If \c n is zero, this function will call the error handler.
966  */
967  template<class vec_t, class vec2_t>
968  double vector_correlation(size_t n, const vec_t &data1,
969  const vec2_t &data2) {
970  size_t i;
971 
972  if (n<1) {
973  O2SCL_ERR2("Cannot compute correlation with no elements",
974  " in vector_correlation().",exc_einval);
975  }
976 
977  double sum_xsq=0.0;
978  double sum_ysq=0.0;
979  double sum_cross=0.0;
980  double ratio;
981  double delta_x, delta_y;
982  double mean_x, mean_y;
983  double r;
984 
985  /*
986  * Compute:
987  * sum_xsq = Sum [ (x_i - mu_x)^2 ],
988  * sum_ysq = Sum [ (y_i - mu_y)^2 ] and
989  * sum_cross = Sum [ (x_i - mu_x) * (y_i - mu_y) ]
990  * using the above relation from Welford's paper
991  */
992 
993  mean_x=data1[0];
994  mean_y=data2[0];
995 
996  for (i=1; i < n; ++i) {
997  ratio=i / (i + 1.0);
998  delta_x=data1[i] - mean_x;
999  delta_y=data2[i] - mean_y;
1000  sum_xsq += delta_x * delta_x * ratio;
1001  sum_ysq += delta_y * delta_y * ratio;
1002  sum_cross += delta_x * delta_y * ratio;
1003  mean_x += delta_x / (i + 1.0);
1004  mean_y += delta_y / (i + 1.0);
1005  }
1006 
1007  r=sum_cross / (std::sqrt(sum_xsq) * std::sqrt(sum_ysq));
1008 
1009  return r;
1010  }
1011 
1012  /** \brief Pearson's correlation
1013 
1014  This function computes the Pearson correlation coefficient
1015  between \c data1 and \c data2 .
1016 
1017  This function produces the same
1018  results as <tt>gsl_stats_correlation()</tt>.
1019 
1020  \comment
1021  r = cov(x, y) / (\Hat\sigma_x \Hat\sigma_y)
1022  = {1/(n-1) \sum (x_i - \Hat x) (y_i - \Hat y)
1023  \over
1024  \sqrt{1/(n-1) \sum (x_i - \Hat x)^2} \sqrt{1/(n-1)
1025  \sum (y_i - \Hat y)^2}
1026  }
1027  \endcomment
1028 
1029  If \c n is zero, this function will call the error handler.
1030  */
1031  template<class vec_t, class vec2_t>
1032  double vector_correlation(const vec_t &data1,
1033  const vec2_t &data2) {
1034  return vector_correlation(data1.size(),data1,data2);
1035  }
1036 
1037  /** \brief The pooled variance of two vectors
1038 
1039  This function computes
1040  \f[
1041  s_{p}^2 = \frac{(n_1-1)s_1^2+(n_2-1)s_2^2}{n_1+n_2-2}
1042  \f]
1043  where \f$ n_i \f$ is the number of elements in vector \f$ i \f$
1044  and \f$ s_i^2 \f$ is the variance of vector \f$ i \f$.
1045 
1046  From http://en.wikipedia.org/wiki/Pooled_variance, "Under the
1047  assumption of equal population variances, the pooled sample
1048  variance provides a higher precision estimate of variance than
1049  the individual sample variances."
1050 
1051  This function produces the same
1052  results as <tt>gsl_stats_pvariance()</tt>.
1053  */
1054  template<class vec_t, class vec2_t>
1055  double vector_pvariance(size_t n1, const vec_t &data1,
1056  size_t n2, const vec2_t &data2) {
1057  double var1=vector_variance<vec_t>(n1,data1);
1058  double var2=vector_variance<vec2_t>(n2,data2);
1059  return (((n1-1)*var1)+((n2-1)*var2))/(n1+n2-2);
1060  }
1061 
1062  /** \brief The pooled variance of two vectors
1063 
1064  This function computes
1065  \f[
1066  s_{p}^2 = \frac{(n_1-1)s_1^2+(n_2-1)s_2^2}{n_1+n_2-2}
1067  \f]
1068  where \f$ n_i \f$ is the number of elements in vector \f$ i \f$
1069  and \f$ s_i^2 \f$ is the variance of vector \f$ i \f$.
1070 
1071  From http://en.wikipedia.org/wiki/Pooled_variance, "Under the
1072  assumption of equal population variances, the pooled sample
1073  variance provides a higher precision estimate of variance than
1074  the individual sample variances."
1075 
1076  This function produces the same
1077  results as <tt>gsl_stats_pvariance()</tt>.
1078  */
1079  template<class vec_t, class vec2_t>
1080  double vector_pvariance(const vec_t &data1,
1081  const vec2_t &data2) {
1082  return vector_pvariance(data1.size(),data1,data2.size(),data2);
1083  }
1084 
1085  /** \brief Quantile from sorted data (ascending only)
1086 
1087  This function returns the quantile \c f of data which
1088  has already been sorted in ascending order. The quantile,
1089  \f$ q \f$ , is
1090  found by interpolation using
1091  \f[
1092  q = \left(1-\delta\right) x_i \delta x_{i+1}
1093  \f]
1094  where \f$ i = \mathrm{floor}[ (n-1)f ] \f$ and
1095  \f$ \delta = (n-1)f -i \f$ .
1096 
1097  This function produces the same
1098  results as <tt>gsl_stats_quantile_from_sorted_data()</tt>.
1099 
1100  No checks are made to ensure the data is sorted, or to ensure
1101  that \f$ 0 \leq 0 \leq 1 \f$. If \c n is zero, this function
1102  will return zero without calling the error handler.
1103  */
1104  template<class vec_t>
1105  double vector_quantile_sorted(size_t n, const vec_t &data,
1106  const double f) {
1107 
1108  double index=f*(n-1);
1109  size_t lhs=((size_t)index);
1110  double delta=index-lhs;
1111  if (n==0) return 0.0;
1112  if (lhs==n-1) return data[lhs];
1113  return (1-delta)*data[lhs]+delta*data[lhs+1];
1114  }
1115 
1116  /** \brief Quantile from sorted data (ascending only)
1117 
1118  This function returns the quantile \c f of data which
1119  has already been sorted in ascending order. The quantile,
1120  \f$ q \f$ , is
1121  found by interpolation using
1122  \f[
1123  q = \left(1-\delta\right) x_i \delta x_{i+1}
1124  \f]
1125  where \f$ i = \mathrm{floor}[ (n-1)f ] \f$ and
1126  \f$ \delta = (n-1)f -i \f$ .
1127 
1128  This function produces the same
1129  results as <tt>gsl_stats_quantile_from_sorted_data()</tt>.
1130 
1131  No checks are made to ensure the data is sorted, or to ensure
1132  that \f$ 0 \leq 0 \leq 1 \f$. If \c n is zero, this function
1133  will return zero without calling the error handler.
1134  */
1135  template<class vec_t>
1136  double vector_quantile_sorted(const vec_t &data, const double f) {
1137  return vector_quantile_sorted<vec_t>(data.size(),data,f);
1138  }
1139 
1140  /** \brief Return the median of sorted (ascending or descending) data
1141 
1142  This function returns the median of sorted data (either
1143  ascending or descending), assuming the data has already been
1144  sorted. When the data set has an odd number of elements, the
1145  median is the value of the element at index \f$ (n-1)/2 \f$,
1146  otherwise, the median is taken to be the average of the elements
1147  at indices \f$ (n-1)/2 \f$ and \f$ n/2 \f$ .
1148 
1149  This function produces the same
1150  results as <tt>gsl_stats_median_from_sorted_data()</tt>.
1151 
1152  No checks are made to ensure the data is sorted. If \c n is
1153  zero, this function will return zero without calling the error
1154  handler.
1155  */
1156  template<class vec_t>
1157  double vector_median_sorted(size_t n, const vec_t &data) {
1158 
1159  if (n==0) return 0.0;
1160 
1161  size_t lhs=(n-1)/2;
1162  size_t rhs=n/2;
1163 
1164  if (lhs==rhs) return data[lhs];
1165 
1166  return (data[lhs]+data[rhs])/2.0;
1167  }
1168 
1169  /** \brief Return the median of sorted (ascending or descending) data
1170 
1171  This function returns the median of sorted data (either
1172  ascending or descending), assuming the data has already been
1173  sorted. When the data set has an odd number of elements, the
1174  median is the value of the element at index \f$ (n-1)/2 \f$,
1175  otherwise, the median is taken to be the average of the elements
1176  at indices \f$ (n-1)/2 \f$ and \f$ n/2 \f$ .
1177 
1178  This function produces the same
1179  results as <tt>gsl_stats_median_from_sorted_data()</tt>.
1180 
1181  No checks are made to ensure the data is sorted. If \c n is
1182  zero, this function will return zero without calling the error
1183  handler.
1184  */
1185  template<class vec_t>
1186  double vector_median_sorted(const vec_t &data) {
1187  return vector_median_sorted<vec_t>(data.size(),data);
1188  }
1189 
1190  /** \brief Compute the chi-squared statistic
1191 
1192  This function computes
1193  \f[
1194  \sum_i \left( \frac{\mathrm{obs}_i - \mathrm{exp}_i}
1195  {\mathrm{err}_i}\right)^2
1196  \f]
1197  where \f$ \mathrm{obs} \f$ are the observed values,
1198  \f$ \mathrm{exp} \f$ are the expected values, and
1199  \f$ \mathrm{err} \f$ are the errors.
1200  */
1201  template<class vec_t, class vec2_t, class vec3_t>
1202  double vector_chi_squared(size_t n, const vec_t &obs, const vec2_t &exp,
1203  const vec3_t &err) {
1204  double chi2=0.0;
1205  for(size_t i=0;i<n;i++) {
1206  chi2+=pow((obs[i]-exp[i])/err[i],2.0);
1207  }
1208  return chi2;
1209  }
1210 
1211  /** \brief Compute the chi-squared statistic
1212 
1213  This function computes
1214  \f[
1215  \sum_i \left( \frac{\mathrm{obs}_i - \mathrm{exp}_i}
1216  {\mathrm{err}_i}\right)^2
1217  \f]
1218  where \f$ \mathrm{obs} \f$ are the observed values,
1219  \f$ \mathrm{exp} \f$ are the expected values, and
1220  \f$ \mathrm{err} \f$ are the errors.
1221  */
1222  template<class vec_t, class vec2_t, class vec3_t>
1223  double vector_chi_squared(const vec_t &obs, const vec2_t &exp,
1224  const vec3_t &err) {
1225  return vector_chi_squared<vec_t,vec2_t,vec3_t>(obs.size(),obs,exp,err);
1226  }
1227 
1228  /** \brief Optimal bin size using Scott's method for the
1229  first <tt>n</tt> elements
1230  */
1231  template<class vec_t> double vector_bin_size_scott
1232  (size_t n, const vec_t &v) {
1233  if (n<=1) return 0.0;
1234  double ret=3.5*vector_stddev(n,v)/cbrt(((double)n));
1235  return ret;
1236  }
1237 
1238  /** \brief Optimal bin size using Scott's method
1239 
1240  This function computes the optimal bin size \f$ \Delta_b \f$ of
1241  a histogram using the expression
1242  \f[
1243  \Delta_b = \frac{3.5 \sigma}{n^{1/3}}
1244  \f]
1245 
1246  From \ref Scott79 .
1247 
1248  \note If <tt>n</tt> is less than or equal to 1, this
1249  function returns 0.0 without calling the error handler.
1250  */
1251  template<class vec_t> double vector_bin_size_scott
1252  (const vec_t &v) {
1253  return vector_bin_size_scott(v.size(),v);
1254  }
1255 
1256  /** \brief Obtain a quantile from a sorted vector
1257 
1258  This is a generic version of
1259  <tt>gsl_stats_quantile_from_sorted_data()</tt>.
1260 
1261  If <tt>f</tt> is less than 0 or greater than 1, the error
1262  handler is called. If <tt>n</tt> is zero, this function returns
1263  zero without calling the error handler.
1264  */
1265  template<class vec_t> double vector_sorted_quantile
1266  (size_t n, const vec_t &v, double f) {
1267 
1268  if (f<0.0 || f>1.0) {
1269  O2SCL_ERR("Invalid fraction for vector_sorted_quantile",
1271  }
1272 
1273  double index=f*(n-1);
1274  size_t lhs=(int)index;
1275  double delta=index-lhs;
1276  double result;
1277 
1278  if (n == 0) {
1279  return 0.0;
1280  }
1281 
1282  if (lhs == n-1) {
1283  result=v[lhs];
1284  } else {
1285  result=(1-delta)*v[lhs]+delta*v[(lhs+1)];
1286  }
1287 
1288  return result;
1289  }
1290 
1291  /** \brief Optimal bin size using the Freedman-Diaconis rule
1292  for the first <tt>n</tt> elements
1293  */
1294  template<class vec_t> double vector_bin_size_freedman
1295  (size_t n, vec_t &v) {
1296  vector_sort<vec_t,double>(n,v);
1297  double ret=2.0*(vector_sorted_quantile(n,v,0.75)-
1298  vector_sorted_quantile(n,v,0.25))/cbrt(((double)n));
1299  return ret;
1300  }
1301 
1302  /** \brief Optimal bin size using the Freedman-Diaconis rule
1303 
1304  This function computes the optimal bin size \f$ \Delta_b \f$ of
1305  a histogram using the expression
1306  \f[
1307  \Delta_b = \frac{2\left(q_{0.75}-q_{0.25}\right)}{n^{1/3}}
1308  \f]
1309  where \f$ q_{i} \f$ is the \f$ i \f$ quantile
1310  of the data (note this is quantile not quartile).
1311  This function sorts the vector in order to obtain
1312  the result.
1313 
1314  From \ref Freedman81 .
1315 
1316  \note If <tt>n</tt> is less than or equal to 1, this
1317  function returns 0.0 without calling the error handler.
1318  */
1319  template<class vec_t> double vector_bin_size_freedman
1320  (vec_t &v) {
1321  return vector_bin_size_freedman(v.size(),v);
1322  }
1323  //@}
1324 
1325  /// \name Weighted vector functions
1326  //@{
1327  /** \brief Compute the mean of weighted data
1328 
1329  This function computes
1330  \f[
1331  \left( \sum_i w_i x_i \right) \left( \sum_i w_i \right)^{-1}
1332  \f]
1333 
1334  This function produces the same results
1335  as <tt>gsl_stats_wmean()</tt>.
1336 
1337  \comment
1338  M(n) = M(n-1) + (data[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
1339  W(n) = W(n-1) + w(n)
1340  \endcomment
1341  */
1342  template<class vec_t, class vec2_t>
1343  double wvector_mean(size_t n, const vec_t &data, const vec2_t &weights) {
1344 
1345  long double wmean=0.0;
1346  long double W=0.0;
1347  for(size_t i=0;i<n;i++) {
1348  double wi=weights[i];
1349  if (wi>0.0) {
1350  W+=wi;
1351  wmean+=(data[i]-wmean)*(wi/W);
1352  }
1353  }
1354 
1355  return wmean;
1356  }
1357 
1358  /** \brief Compute the mean of weighted data
1359 
1360  This function computes
1361  \f[
1362  \left( \sum_i w_i x_i \right) \left( \sum_i w_i \right)^{-1}
1363  \f]
1364 
1365  This function produces the same results
1366  as <tt>gsl_stats_wmean()</tt>.
1367 
1368  \comment
1369  M(n) = M(n-1) + (data[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
1370  W(n) = W(n-1) + w(n)
1371  \endcomment
1372  */
1373  template<class vec_t, class vec2_t>
1374  double wvector_mean(const vec_t &data, const vec2_t &weights) {
1375  return wvector_mean<vec_t,vec2_t>(data.size(),data,weights);
1376  }
1377 
1378  /** \brief Compute a normalization factor for weighted data
1379 
1380  This function is used internally in \ref wvector_variance(size_t
1381  n, vec_t &data, const vec2_t &weights, double wmean) and \ref
1382  wvector_stddev(size_t n, vec_t &data, const vec2_t &weights, double
1383  wmean) .
1384  */
1385  template<class vec_t> double wvector_factor(size_t n, const vec_t &weights) {
1386 
1387  long double a=0.0;
1388  long double b=0.0;
1389  long double factor;
1390  for(size_t i=0;i<n;i++) {
1391  double wi=weights[i];
1392  if (wi>0.0) {
1393  a+=wi;
1394  b+=wi*wi;
1395  }
1396  }
1397  factor=a*a/(a*a-b);
1398  return factor;
1399  }
1400 
1401  /** \brief Compute a normalization factor for weighted data
1402 
1403  This function is used internally in \ref wvector_variance(size_t
1404  n, vec_t &data, const vec2_t &weights, double wmean) and \ref
1405  wvector_stddev(size_t n, vec_t &data, const vec2_t &weights, double
1406  wmean) .
1407  */
1408  template<class vec_t> double wvector_factor(const vec_t &weights) {
1409  return wvector_factor<vec_t>(weights.size(),weights);
1410  }
1411 
1412  /** \brief Compute the variance of a weighted vector with a mean
1413  known in advance
1414 
1415  This function computes
1416  \f[
1417  \left[ \sum_i w_i \left(x_i-\mu\right)^2 \right]
1418  \left[ \sum_i w_i \right]^{-1}
1419  \f]
1420 
1421  This function produces the same results
1422  as <tt>gsl_stats_wvariance_with_fixed_mean()</tt>.
1423 
1424  */
1425  template<class vec_t, class vec2_t>
1426  double wvector_variance_fmean(size_t n, const vec_t &data,
1427  const vec2_t &weights, double wmean) {
1428  long double wvariance=0.0;
1429  long double W=0.0;
1430  for(size_t i=0;i<n;i++) {
1431  double wi=weights[i];
1432  if (wi>0.0) {
1433  const long double delta=data[i]-wmean;
1434  W+=wi;
1435  wvariance+=(delta*delta-wvariance)*(wi/W);
1436  }
1437  }
1438 
1439  return wvariance;
1440  }
1441 
1442  /** \brief Compute the variance of a weighted vector with a mean
1443  known in advance
1444 
1445  This function computes
1446  \f[
1447  \left[ \sum_i w_i \left(x_i-\mu\right)^2 \right]
1448  \left[ \sum_i w_i \right]^{-1}
1449  \f]
1450 
1451  This function produces the same results
1452  as <tt>gsl_stats_wvariance_with_fixed_mean()</tt>.
1453 
1454  */
1455  template<class vec_t, class vec2_t>
1456  double wvector_variance_fmean(const vec_t &data,
1457  const vec2_t &weights, double wmean) {
1458  return wvector_variance_fmean(data.size(),data,weights,wmean);
1459  }
1460 
1461  /** \brief Compute the variance of a weighted vector with
1462  specified mean
1463 
1464  This function produces the same results
1465  as <tt>gsl_stats_wvariance_m()</tt>.
1466  */
1467  template<class vec_t, class vec2_t>
1468  double wvector_variance(size_t n, const vec_t &data,
1469  const vec2_t &weights, double wmean) {
1470 
1471  const double variance=wvector_variance_fmean
1472  (n,data,weights,wmean);
1473  const double scale=wvector_factor(n,weights);
1474  const double wvar=scale*variance;
1475  return wvar;
1476  }
1477 
1478  /** \brief Compute the variance of a weighted vector with
1479  specified mean
1480 
1481  This function produces the same results
1482  as <tt>gsl_stats_wvariance_m()</tt>.
1483  */
1484  template<class vec_t, class vec2_t>
1485  double wvector_variance(const vec_t &data,
1486  const vec2_t &weights, double wmean) {
1487  return wvector_variance<vec_t,vec2_t>(data.size(),data,weights,wmean);
1488  }
1489 
1490  /** \brief Compute the variance of a weighted vector where mean
1491  is computed automatically
1492 
1493  This function produces the same results
1494  as <tt>gsl_stats_wvariance()</tt>.
1495  */
1496  template<class vec_t, class vec2_t>
1497  double wvector_variance(size_t n, const vec_t &data,
1498  const vec2_t &weights) {
1499 
1500  double wmean=wvector_mean(n,data,weights);
1501  return wvector_variance<vec_t,vec2_t>(n,data,weights,wmean);
1502  }
1503 
1504  /** \brief Compute the variance of a weighted vector where mean
1505  is computed automatically
1506 
1507  This function produces the same results
1508  as <tt>gsl_stats_wvariance()</tt>.
1509  */
1510  template<class vec_t, class vec2_t>
1511  double wvector_variance(const vec_t &data, const vec2_t &weights) {
1512  return wvector_variance(data.size(),data,weights);
1513  }
1514 
1515  /** \brief The weighted covariance of two vectors
1516 
1517  \note Experimental
1518  */
1519  template<class vec_t, class vec2_t, class vec3_t>
1520  double wvector_covariance(size_t n, const vec_t &data1,
1521  const vec2_t &data2,
1522  const vec3_t &weights) {
1523  double mean1=wvector_mean(n,data1,weights);
1524  double mean2=wvector_mean(n,data2,weights);
1525  double covar=0.0;
1526  double W=0.0;
1527  for(size_t i=0;i<n;i++) {
1528  double wi=weights[i];
1529  if (wi>0.0) {
1530  W+=wi;
1531  double delta1=(data1[i]-mean1);
1532  double delta2=(data2[i]-mean2);
1533  covar+=(wi/W)*(delta1*delta2-covar);
1534  }
1535  }
1536  double scale=wvector_factor(n,weights);
1537  return covar*scale;
1538  }
1539 
1540  /** \brief The weighted covariance of two vectors
1541 
1542  \note Experimental
1543  */
1544  template<class vec_t, class vec2_t, class vec3_t>
1545  double wvector_covariance(const vec_t &data1, const vec2_t &data2,
1546  const vec3_t &weights) {
1547  return wvector_covariance<vec_t,vec2_t,vec3_t>
1548  (data1.size(),data1,data2,weights);
1549  }
1550 
1551  /** \brief Compute the standard deviation of a weighted vector
1552  with a mean known in advance
1553 
1554  This function produces the same results
1555  as <tt>gsl_stats_wsd_with_fixed_mean()</tt>.
1556  */
1557  template<class vec_t, class vec2_t>
1558  double wvector_stddev_fmean(size_t n, const vec_t &data,
1559  const vec2_t &weights, double wmean) {
1560  return sqrt(wvector_variance_fmean(n,data,weights,wmean));
1561  }
1562 
1563  /** \brief Compute the standard deviation of a weighted vector
1564  with a mean known in advance
1565 
1566  This function produces the same results
1567  as <tt>gsl_stats_wsd_with_fixed_mean()</tt>.
1568  */
1569  template<class vec_t, class vec2_t>
1570  double wvector_stddev_fmean(const vec_t &data,
1571  const vec2_t &weights, double wmean) {
1572  return wvector_stddev_fmean<vec_t,vec2_t>
1573  (data.size(),data,weights,wmean);
1574  }
1575 
1576  /** \brief Compute the standard deviation of a weighted vector where mean
1577  is computed automatically
1578 
1579  This function produces the same results
1580  as <tt>gsl_stats_wsd()</tt>.
1581  */
1582  template<class vec_t, class vec2_t>
1583  double wvector_stddev(size_t n, const vec_t &data,
1584  const vec2_t &weights) {
1585  double wmean=wvector_mean(n,data,weights);
1586  return sqrt(wvector_variance(n,data,weights,wmean));
1587  }
1588 
1589  /** \brief Compute the standard deviation of a weighted vector where mean
1590  is computed automatically
1591 
1592  This function produces the same results
1593  as <tt>gsl_stats_wsd()</tt>.
1594  */
1595  template<class vec_t, class vec2_t>
1596  double wvector_stddev(const vec_t &data, const vec2_t &weights) {
1597  return wvector_stddev(data.size(),data,weights);
1598  }
1599 
1600  /** \brief Compute the standard deviation of a weighted vector with
1601  specified mean
1602 
1603  This function produces the same results
1604  as <tt>gsl_stats_wsd_m()</tt>.
1605  */
1606  template<class vec_t, class vec2_t>
1607  double wvector_stddev(size_t n, const vec_t &data,
1608  const vec2_t &weights, double wmean) {
1609  const double variance=wvector_variance_fmean
1610  (n,data,weights,wmean);
1611  const double scale=wvector_factor(n,weights);
1612  double wvar=scale*variance;
1613  return sqrt(wvar);
1614  }
1615 
1616  /** \brief Compute the standard deviation of a weighted vector with
1617  specified mean
1618 
1619  This function produces the same results
1620  as <tt>gsl_stats_wsd_m()</tt>.
1621  */
1622  template<class vec_t, class vec2_t>
1623  double wvector_stddev(const vec_t &data,
1624  const vec2_t &weights, double wmean) {
1625  return wvector_stddev<vec_t,vec2_t>(data.size(),data,weights,wmean);
1626  }
1627 
1628  /** \brief Compute the weighted sum of squares of data about the
1629  specified weighted mean
1630 
1631  This function produces the same results
1632  as <tt>gsl_stats_wtss_m()</tt>.
1633  */
1634  template<class vec_t, class vec2_t>
1635  double wvector_sumsq(size_t n, const vec_t &data,
1636  const vec2_t &weights, double wmean) {
1637  long double wtss=0.0;
1638  for(size_t i=0;i<n;i++) {
1639  double wi=weights[i];
1640  if (wi>0.0) {
1641  const long double delta=data[i]-wmean;
1642  wtss+=wi*delta*delta;
1643  }
1644  }
1645 
1646  return wtss;
1647  }
1648 
1649  /** \brief Compute the weighted sum of squares of data about the
1650  specified weighted mean
1651 
1652  This function produces the same results
1653  as <tt>gsl_stats_wtss_m()</tt>.
1654  */
1655  template<class vec_t, class vec2_t>
1656  double wvector_sumsq(const vec_t &data,
1657  const vec2_t &weights, double wmean) {
1658  return wvector_sumsq<vec_t,vec2_t>(data.size(),data,weights,wmean);
1659  }
1660 
1661  /** \brief Compute the weighted sum of squares of data about the
1662  weighted mean
1663 
1664  This function produces the same results
1665  as <tt>gsl_stats_wtss()</tt>.
1666  */
1667  template<class vec_t, class vec2_t>
1668  double wvector_sumsq(size_t n, const vec_t &data,
1669  const vec2_t &weights) {
1670 
1671  double wmean=wvector_mean(n,data,weights);
1672  return wvector_sumsq(n,data,weights,wmean);
1673  }
1674 
1675  /** \brief Compute the weighted sum of squares of data about the
1676  weighted mean
1677 
1678  This function produces the same results
1679  as <tt>gsl_stats_wtss()</tt>.
1680  */
1681  template<class vec_t, class vec2_t>
1682  double wvector_sumsq(const vec_t &data, const vec2_t &weights) {
1683  return wvector_sumsq<vec_t,vec2_t>(data.size(),data,weights);
1684  }
1685 
1686  /** \brief Compute the absolute deviation of data about a specified mean
1687 
1688  This function produces the same results
1689  as <tt>gsl_stats_wabsdev_m()</tt>.
1690  */
1691  template<class vec_t, class vec2_t>
1692  double wvector_absdev(size_t n, const vec_t &data, const vec2_t &weights,
1693  double wmean) {
1694  long double wabsdev=0.0;
1695  long double W=0.0;
1696  for(size_t i=0;i<n;i++) {
1697  double wi=weights[i];
1698  if (wi>0.0) {
1699  const long double delta=fabs(data[i]-wmean);
1700  W+=wi;
1701  wabsdev+=(delta-wabsdev)*(wi/W);
1702  }
1703  }
1704  return wabsdev;
1705  }
1706 
1707  /** \brief Compute the absolute deviation of data about a specified mean
1708 
1709  This function produces the same results
1710  as <tt>gsl_stats_wabsdev_m()</tt>.
1711  */
1712  template<class vec_t, class vec2_t>
1713  double wvector_absdev(const vec_t &data, const vec2_t &weights,
1714  double wmean) {
1715  return wvector_absdev<vec_t,vec2_t>(data.size(),data,weights,wmean);
1716  }
1717 
1718  /** \brief Compute the absolute deviation of data about a specified mean
1719 
1720  This function produces the same results
1721  as <tt>gsl_stats_wabsdev()</tt>.
1722  */
1723  template<class vec_t, class vec2_t>
1724  double wvector_absdev(size_t n, const vec_t &data,
1725  const vec2_t &weights) {
1726 
1727  double wmean=wvector_mean(n,data,weights);
1728  return wvector_absdev(n,data,weights,wmean);
1729  }
1730 
1731  /** \brief Compute the absolute deviation of data about a specified mean
1732 
1733  This function produces the same results
1734  as <tt>gsl_stats_wabsdev()</tt>.
1735  */
1736  template<class vec_t, class vec2_t>
1737  double wvector_absdev(const vec_t &data, const vec2_t &weights) {
1738  return wvector_absdev<vec_t,vec2_t>(data.size(),data,weights);
1739  }
1740 
1741  /** \brief Compute the skewness of data with specified mean
1742  and standard deviation
1743 
1744  This function produces the same results
1745  as <tt>gsl_stats_wskew_m_sd()</tt>.
1746  */
1747  template<class vec_t, class vec2_t>
1748  double wvector_skew(size_t n, const vec_t &data, const vec2_t &weights,
1749  double wmean, double wsd) {
1750  long double wskew=0.0;
1751  long double W=0.0;
1752  for(size_t i=0;i<n;i++) {
1753  double wi=weights[i];
1754  if (wi>0.0) {
1755  const long double x=(data[i]-wmean)/wsd;
1756  W+=wi;
1757  wskew+=(x*x*x-wskew)*(wi/W);
1758  }
1759  }
1760  return wskew;
1761  }
1762 
1763  /** \brief Compute the skewness of data with specified mean
1764  and standard deviation
1765 
1766  This function produces the same results
1767  as <tt>gsl_stats_wskew_m_sd()</tt>.
1768  */
1769  template<class vec_t, class vec2_t>
1770  double wvector_skew(const vec_t &data, const vec2_t &weights,
1771  double wmean, double wsd) {
1772  return wvector_skew<vec_t,vec2_t>(data.size(),data,weights,wmean,wsd);
1773  }
1774 
1775  /** \brief Compute the skewness of data with specified mean
1776  and standard deviation
1777 
1778  This function produces the same results
1779  as <tt>gsl_stats_wskew()</tt>.
1780  */
1781  template<class vec_t, class vec2_t>
1782  double wvector_skew(size_t n, const vec_t &data, const vec2_t &weights) {
1783  double wmean=wvector_mean(n,data,weights);
1784  double wsd=wvector_stddev(n,data,weights,wmean);
1785  return wvector_skew(n,data,weights,wmean,wsd);
1786  }
1787 
1788  /** \brief Compute the skewness of data with specified mean
1789  and standard deviation
1790 
1791  This function produces the same results
1792  as <tt>gsl_stats_wskew()</tt>.
1793  */
1794  template<class vec_t, class vec2_t>
1795  double wvector_skew(const vec_t &data, const vec2_t &weights) {
1796  return wvector_skew<vec_t,vec2_t>(data.size(),data,weights);
1797  }
1798 
1799  /** \brief Compute the kurtosis of data with specified mean
1800  and standard deviation
1801 
1802  This function produces the same results
1803  as <tt>gsl_stats_wkurtosis_m_sd()</tt>.
1804  */
1805  template<class vec_t, class vec2_t>
1806  double wvector_kurtosis(size_t n, const vec_t &data, const vec2_t &weights,
1807  double wmean, double wsd) {
1808  long double wavg=0.0;
1809  long double W=0.0;
1810  for(size_t i=0;i<n;i++) {
1811  double wi=weights[i];
1812  if (wi>0.0) {
1813  const long double x=(data[i]-wmean)/wsd;
1814  W+=wi;
1815  wavg+=(x*x*x*x-wavg)*(wi/W);
1816  }
1817  }
1818  return wavg-3.0;
1819  }
1820 
1821  /** \brief Compute the kurtosis of data with specified mean
1822  and standard deviation
1823 
1824  This function produces the same results
1825  as <tt>gsl_stats_wkurtosis_m_sd()</tt>.
1826  */
1827  template<class vec_t, class vec2_t>
1828  double wvector_kurtosis(const vec_t &data, const vec2_t &weights,
1829  double wmean, double wsd) {
1830  return wvector_kurtosis<vec_t,vec2_t>
1831  (data.size(),data,weights,wmean,wsd);
1832  }
1833 
1834  /** \brief Compute the kurtosis of data with specified mean
1835  and standard deviation
1836 
1837  This function produces the same results
1838  as <tt>gsl_stats_wkurtosis()</tt>.
1839  */
1840  template<class vec_t, class vec2_t>
1841  double wvector_kurtosis(size_t n, const vec_t &data,
1842  const vec2_t &weights) {
1843  double wmean=wvector_mean(n,data,weights);
1844  double wsd=wvector_stddev(n,data,weights,wmean);
1845  return wvector_kurtosis(n,data,weights,wmean,wsd);
1846  }
1847 
1848  /** \brief Compute the kurtosis of data with specified mean
1849  and standard deviation
1850 
1851  This function produces the same results
1852  as <tt>gsl_stats_wkurtosis()</tt>.
1853  */
1854  template<class vec_t, class vec2_t>
1855  double wvector_kurtosis(const vec_t &data, const vec2_t &weights) {
1856  return wvector_kurtosis<vec_t,vec2_t>(data,weights);
1857  }
1858  //@}
1859 
1860  /** \brief Lag-k autocorrelation for the first
1861  \c n elements with a vector multiplier given
1862  the mean
1863  */
1864  template<class vec_t, class vec2_t>
1865  double vector_lagk_autocorr_mult(size_t n, const vec_t &data,
1866  const vec2_t &mult, size_t k,
1867  double mean) {
1868 
1869  size_t n2=0;
1870  for(size_t i=0;i<n;i++) {
1871  size_t m=((size_t)(mult[i]*(1.0+1.0e-10)));
1872  if (m==0) {
1873  O2SCL_ERR2("Mult vector is zero ",
1874  "in vector_lagk_autocorr_mult().",exc_einval);
1875  }
1876  n2+=m;
1877  }
1878 
1879  if (n2<=k) {
1880  O2SCL_ERR2("Not enough elements ",
1881  "in vector_lagk_autocorr_mult().",exc_einval);
1882  }
1883 
1884  long double q=0.0, v=0.0;
1885  size_t im=0, ix=0, im2=0, ix2=0;
1886  for(size_t i=0;i<k;i++) {
1887  v+=(data[ix]-mean)*(data[ix]-mean)/(i+1);
1888  im++;
1889  if (im>=((size_t)(mult[ix]*(1.0+1.0e-10)))) {
1890  im=0;
1891  ix++;
1892  }
1893  }
1894  for(size_t i=k;i<n2;i++) {
1895  long double delta0=data[ix2]-mean;
1896  long double delta1=data[ix]-mean;
1897  q+=(delta0*delta1-q)/(i+1);
1898  v+=(delta1*delta1-v)/(i+1);
1899  im++;
1900  if (im>=((size_t)(mult[ix]*(1.0+1.0e-10)))) {
1901  im=0;
1902  ix++;
1903  }
1904  im2++;
1905  if (im2>=((size_t)(mult[ix2]*(1.0+1.0e-10)))) {
1906  im2=0;
1907  ix2++;
1908  }
1909  }
1910  return q/v;
1911  }
1912 
1913  /** \brief Lag-k autocorrelation for the first
1914  \c n elements with a vector multiplier
1915  */
1916  template<class vec_t, class vec2_t>
1917  double vector_lagk_autocorr_mult(size_t n, const vec_t &data,
1918  const vec2_t &mult, size_t k) {
1919  double mean=wvector_mean(n,mult,data);
1920  return vector_lagk_autocorr_mult(n,data,mult,k,mean);
1921  }
1922 
1923  /** \brief Lag-k autocorrelation with a vector multiplier
1924  given the mean
1925  */
1926  template<class vec_t, class vec2_t>
1927  double vector_lagk_autocorr_mult(const vec_t &data,
1928  const vec2_t &mult, size_t k,
1929  double mean) {
1930  return vector_lagk_autocorr_mult(data.size(),data,mult,k,mean);
1931  }
1932 
1933  /** \brief Lag-k autocorrelation with a vector multiplier
1934  */
1935  template<class vec_t, class vec2_t>
1936  double vector_lagk_autocorr_mult(const vec_t &data,
1937  const vec2_t &mult, size_t k) {
1938  return vector_lagk_autocorr_mult(data.size(),data,mult,k);
1939  }
1940 
1941  /** \brief Construct an autocorrelation vector using a multiplier
1942  using the first \c n2 elements of vectors \c data and \c mult
1943  */
1944  template<class vec_t, class vec2_t, class resize_vec_t>
1946  (size_t n2, const vec_t &data, const vec2_t &mult, resize_vec_t &ac_vec) {
1947 
1948  size_t n=0;
1949  for(size_t i=0;i<n2;i++) {
1950  n+=((size_t)(mult[i]*(1.0+1.0e-10)));
1951  }
1952 
1953  size_t kmax=n/2;
1954  double mean=wvector_mean(data,mult);
1955  ac_vec.resize(kmax);
1956  ac_vec[0]=1.0;
1957  for(size_t k=1;k<kmax;k++) {
1958  ac_vec[k]=vector_lagk_autocorr_mult(n2,data,mult,k,mean);
1959  }
1960  return;
1961  }
1962 
1963  /** \brief Construct an autocorrelation vector using a multiplier
1964  */
1965  template<class vec_t, class vec2_t, class resize_vec_t>
1967  (const vec_t &data, const vec2_t &mult, resize_vec_t &ac_vec) {
1968 
1969  vector_autocorr_vector_mult(data.size(),data,mult,ac_vec);
1970 
1971  return;
1972  }
1973 
1974 #ifndef DOXYGEN_NO_O2NS
1975 }
1976 #endif
1977 
1978 #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 vector_sorted_quantile(size_t n, const vec_t &v, double f)
Obtain a quantile from a sorted vector.
Definition: vec_stats.h:1266
double wvector_mean(size_t n, const vec_t &data, const vec2_t &weights)
Compute the mean of weighted data.
Definition: vec_stats.h:1343
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:1692
void vector_autocorr_vector(const vec_t &data, resize_vec_t &ac_vec)
Construct an autocorrelation vector.
Definition: vec_stats.h:788
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:1426
double vector_bin_size_scott(size_t n, const vec_t &v)
Optimal bin size using Scott&#39;s method for the first n elements.
Definition: vec_stats.h:1232
double vector_correlation(size_t n, const vec_t &data1, const vec2_t &data2)
Pearson&#39;s correlation.
Definition: vec_stats.h:968
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:1806
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:1946
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:1055
double vector_median_sorted(size_t n, const vec_t &data)
Return the median of sorted (ascending or descending) data.
Definition: vec_stats.h:1157
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:1558
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:1520
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:1583
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:1468
#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:1635
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
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:861
double vector_quantile_sorted(size_t n, const vec_t &data, const double f)
Quantile from sorted data (ascending only)
Definition: vec_stats.h:1105
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:1202
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:1865
double wvector_factor(size_t n, const vec_t &weights)
Compute a normalization factor for weighted data.
Definition: vec_stats.h:1385
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:1295
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:1748
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
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:825

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