yat/statistics/utility.cc

Code
Comments
Other
Rev Date Author Line
115 19 Jul 04 peter 1 // $Id$
115 19 Jul 04 peter 2
675 10 Oct 06 jari 3 /*
2119 12 Dec 09 peter 4   Copyright (C) 2004 Jari Häkkinen, Peter Johansson
831 27 Mar 07 peter 5   Copyright (C) 2005 Peter Johansson
2119 12 Dec 09 peter 6   Copyright (C) 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 7   Copyright (C) 2009, 2011, 2012, 2018, 2021, 2022, 2023 Peter Johansson
115 19 Jul 04 peter 8
1437 25 Aug 08 peter 9   This file is part of the yat library, http://dev.thep.lu.se/yat
298 29 Apr 05 peter 10
675 10 Oct 06 jari 11   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 12   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 13   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 14   License, or (at your option) any later version.
675 10 Oct 06 jari 15
675 10 Oct 06 jari 16   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
675 10 Oct 06 jari 19   General Public License for more details.
675 10 Oct 06 jari 20
675 10 Oct 06 jari 21   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 22   along with yat. If not, see <http://www.gnu.org/licenses/>.
675 10 Oct 06 jari 23 */
675 10 Oct 06 jari 24
2881 18 Nov 12 peter 25 #include <config.h>
2881 18 Nov 12 peter 26
680 11 Oct 06 jari 27 #include "utility.h"
675 10 Oct 06 jari 28
1145 25 Feb 08 peter 29 #include <gsl/gsl_cdf.h>
169 23 Sep 04 peter 30 #include <gsl/gsl_randist.h>
4336 14 Apr 23 peter 31 #include <gsl/gsl_sf_gamma.h>
500 28 Jan 06 peter 32 #include <gsl/gsl_statistics_double.h>
115 19 Jul 04 peter 33
1746 23 Jan 09 peter 34 #include <algorithm>
1145 25 Feb 08 peter 35 #include <cassert>
2641 10 Nov 11 peter 36 #include <cmath>
2055 08 Sep 09 peter 37 #include <limits>
1145 25 Feb 08 peter 38
115 19 Jul 04 peter 39 namespace theplu {
680 11 Oct 06 jari 40 namespace yat {
2641 10 Nov 11 peter 41 namespace statistics {
115 19 Jul 04 peter 42
2641 10 Nov 11 peter 43   double cdf_hypergeometric_P(unsigned int k, unsigned int n1,
1271 09 Apr 08 peter 44                               unsigned int n2, unsigned int t)
169 23 Sep 04 peter 45   {
1287 24 Apr 08 peter 46     return gsl_cdf_hypergeometric_P(k, n1, n2, t);
169 23 Sep 04 peter 47   }
169 23 Sep 04 peter 48
1145 25 Feb 08 peter 49
4125 14 Jan 22 peter 50   utility::Matrix correlation(const utility::MatrixBase& X)
3745 27 Jul 18 peter 51   {
4053 26 Mar 21 peter 52     size_t cols = X.columns();
4053 26 Mar 21 peter 53     size_t rows = X.rows();
4053 26 Mar 21 peter 54     utility::Vector m(cols);
4053 26 Mar 21 peter 55     utility::Vector x2(cols);
4053 26 Mar 21 peter 56     utility::Vector stddev(cols);
4053 26 Mar 21 peter 57     for (size_t i=0; i<cols; ++i) {
4053 26 Mar 21 peter 58       utility::VectorConstView vec = X.column_const_view(i);
4053 26 Mar 21 peter 59       m(i) = sum(vec);
4053 26 Mar 21 peter 60       x2(i) = vec * vec;
4053 26 Mar 21 peter 61       // scaled standard deviation
4053 26 Mar 21 peter 62       stddev(i) = std::sqrt(x2(i) - m(i) * m(i) / rows);
4053 26 Mar 21 peter 63     }
3745 27 Jul 18 peter 64
4053 26 Mar 21 peter 65     utility::Matrix corr(cols, cols, 1.0);
4053 26 Mar 21 peter 66     for (size_t i=0; i<cols; ++i)
4053 26 Mar 21 peter 67       for (size_t j=i+1; j<cols; ++j) {
4053 26 Mar 21 peter 68         corr(i,j) =
4053 26 Mar 21 peter 69           X.column_const_view(i) * X.column_const_view(j) - m(i)*m(j)/rows;
4053 26 Mar 21 peter 70         corr(i, j) /= stddev(i) * stddev(j);
4053 26 Mar 21 peter 71         // symmetry
4053 26 Mar 21 peter 72         corr(j, i) = corr(i,j);
3745 27 Jul 18 peter 73       }
3745 27 Jul 18 peter 74
3745 27 Jul 18 peter 75     return corr;
3745 27 Jul 18 peter 76   }
3745 27 Jul 18 peter 77
3745 27 Jul 18 peter 78
1271 09 Apr 08 peter 79   double pearson_p_value(double r, unsigned int n)
1145 25 Feb 08 peter 80   {
1145 25 Feb 08 peter 81     assert(n>=2);
1271 09 Apr 08 peter 82     if (n<=2)
1145 25 Feb 08 peter 83       return std::numeric_limits<double>::quiet_NaN();
1553 07 Oct 08 peter 84     if (r>=1.0)
1553 07 Oct 08 peter 85       return 0.0;
1553 07 Oct 08 peter 86     if (r<=-1.0)
1553 07 Oct 08 peter 87       return 1.0;
2641 10 Nov 11 peter 88     return gsl_cdf_tdist_Q(r*std::sqrt((n-2)/(1-r*r)), n-2);
1145 25 Feb 08 peter 89   }
1145 25 Feb 08 peter 90
1145 25 Feb 08 peter 91
1025 01 Feb 08 peter 92   double kurtosis(const utility::VectorBase& v)
703 18 Dec 06 jari 93   {
703 18 Dec 06 jari 94     const gsl_vector* gvp=v.gsl_vector_p();
703 18 Dec 06 jari 95     return gsl_stats_kurtosis(gvp->data,gvp->stride,gvp->size);
703 18 Dec 06 jari 96   }
703 18 Dec 06 jari 97
672 07 Oct 06 peter 98
4336 14 Apr 23 peter 99   double negative_hypergeometric_pdf(unsigned k, unsigned n1,
4336 14 Apr 23 peter 100                                      unsigned n2, unsigned t)
4336 14 Apr 23 peter 101   {
4336 14 Apr 23 peter 102     return gsl_sf_choose(k+t+1, k) * gsl_sf_choose(n1+n2-t-k, n1-k)
4336 14 Apr 23 peter 103       / gsl_sf_choose(n1+n2, n1);
4336 14 Apr 23 peter 104   }
4336 14 Apr 23 peter 105
4336 14 Apr 23 peter 106
1025 01 Feb 08 peter 107   double skewness(const utility::VectorBase& v)
703 18 Dec 06 jari 108   {
703 18 Dec 06 jari 109     const gsl_vector* gvp=v.gsl_vector_p();
703 18 Dec 06 jari 110     return gsl_stats_skew(gvp->data,gvp->stride,gvp->size);
703 18 Dec 06 jari 111   }
703 18 Dec 06 jari 112
683 11 Oct 06 jari 113 }}} // of namespace statistics, yat, and theplu

Generated on Sat Oct 14 02:24:08 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)