source: trunk/yat/statistics/utility.cc

Last change on this file was 4089, checked in by Peter, 13 days ago

updating copyright years

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 2.6 KB
Line 
1// $Id: utility.cc 4089 2021-09-07 00:56:40Z peter $
2
3/*
4  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
5  Copyright (C) 2005 Peter Johansson
6  Copyright (C) 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
7  Copyright (C) 2009, 2011, 2012, 2018, 2021 Peter Johansson
8
9  This file is part of the yat library, http://dev.thep.lu.se/yat
10
11  The yat library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU General Public License as
13  published by the Free Software Foundation; either version 3 of the
14  License, or (at your option) any later version.
15
16  The yat library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  General Public License for more details.
20
21  You should have received a copy of the GNU General Public License
22  along with yat. If not, see <http://www.gnu.org/licenses/>.
23*/
24
25#include <config.h>
26
27#include "utility.h"
28
29#include <gsl/gsl_cdf.h>
30#include <gsl/gsl_randist.h>
31#include <gsl/gsl_statistics_double.h>
32
33#include <algorithm>
34#include <cassert>
35#include <cmath>
36#include <limits>
37
38namespace theplu {
39namespace yat {
40namespace statistics {
41
42  double cdf_hypergeometric_P(unsigned int k, unsigned int n1,
43                              unsigned int n2, unsigned int t)
44  {
45    return gsl_cdf_hypergeometric_P(k, n1, n2, t);
46  }
47
48
49  utility::Matrix correlation(const utility::Matrix& X)
50  {
51    size_t cols = X.columns();
52    size_t rows = X.rows();
53    utility::Vector m(cols);
54    utility::Vector x2(cols);
55    utility::Vector stddev(cols);
56    for (size_t i=0; i<cols; ++i) {
57      utility::VectorConstView vec = X.column_const_view(i);
58      m(i) = sum(vec);
59      x2(i) = vec * vec;
60      // scaled standard deviation
61      stddev(i) = std::sqrt(x2(i) - m(i) * m(i) / rows);
62    }
63
64    utility::Matrix corr(cols, cols, 1.0);
65    for (size_t i=0; i<cols; ++i)
66      for (size_t j=i+1; j<cols; ++j) {
67        corr(i,j) =
68          X.column_const_view(i) * X.column_const_view(j) - m(i)*m(j)/rows;
69        corr(i, j) /= stddev(i) * stddev(j);
70        // symmetry
71        corr(j, i) = corr(i,j);
72      }
73
74    return corr;
75  }
76
77
78  double pearson_p_value(double r, unsigned int n)
79  {
80    assert(n>=2);
81    if (n<=2)
82      return std::numeric_limits<double>::quiet_NaN();
83    if (r>=1.0)
84      return 0.0;
85    if (r<=-1.0)
86      return 1.0;
87    return gsl_cdf_tdist_Q(r*std::sqrt((n-2)/(1-r*r)), n-2);
88  }
89
90
91  double kurtosis(const utility::VectorBase& v)
92  {
93    const gsl_vector* gvp=v.gsl_vector_p();
94    return gsl_stats_kurtosis(gvp->data,gvp->stride,gvp->size);
95  }
96
97
98  double skewness(const utility::VectorBase& v)
99  {
100    const gsl_vector* gvp=v.gsl_vector_p();
101    return gsl_stats_skew(gvp->data,gvp->stride,gvp->size);
102  }
103
104}}} // of namespace statistics, yat, and theplu
Note: See TracBrowser for help on using the repository browser.