Changeset 3745


Ignore:
Timestamp:
Jul 27, 2018, 2:27:02 AM (4 years ago)
Author:
Peter
Message:

add function calculating correlation matrix

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/test/statistics.cc

    r3550 r3745  
    5050void test_benjamini_hochberg(test::Suite&);
    5151void test_benjamini_hochberg_unsorted(test::Suite&);
     52void test_correlation(test::Suite& suite);
    5253void test_entropy(test::Suite&);
    5354void test_mad(test::Suite&);
     
    142143  test_median_empty(suite);
    143144  test_mutual_information(suite);
     145  test_correlation(suite);
    144146  return suite.return_value();
     147}
     148
     149
     150void test_correlation(test::Suite& suite)
     151{
     152  utility::Matrix x(10,100);
     153  utility::Matrix C = statistics::correlation(x);
     154  if (!suite.equal(C(0,0), 1.0)) {
     155    suite.add(false);
     156    suite.err() << "correlation element(0, 0) not unity\n";
     157  }
    145158}
    146159
  • trunk/yat/statistics/utility.cc

    r2919 r3745  
    2727#include "utility.h"
    2828
     29#include "yat/normalizer/Centralizer.h"
     30#include "yat/normalizer/ColumnNormalizer.h"
     31
    2932#include <gsl/gsl_cdf.h>
    3033#include <gsl/gsl_randist.h>
     
    4447  {
    4548    return gsl_cdf_hypergeometric_P(k, n1, n2, t);
     49  }
     50
     51
     52  utility::Matrix correlation(utility::Matrix x)
     53  {
     54    using utility::Matrix;
     55    size_t n = x.columns();
     56    // centralise
     57    normalizer::ColumnNormalizer<normalizer::Centralizer<> > normalizer;
     58    normalizer(x, x);
     59
     60    Matrix cov(n, n);
     61    for (size_t i=0; i<n; ++i) {
     62      for (size_t j=i; j<n; ++j) {
     63        cov(i, j) = x.column_const_view(i) * x.column_const_view(j);
     64        cov(j, i) = cov(i, j);
     65      }
     66    }
     67
     68    utility::Vector stddev(n);
     69    for (size_t i=0; i<n; ++i)
     70      stddev(i) = std::sqrt(cov(i, i));
     71
     72    Matrix corr(cov);
     73    for (size_t i=0; i<n; ++i) {
     74      corr(i, i) = 1.0;
     75      for (size_t j=0; j<i; ++j) {
     76        corr(i, j) = cov(i, j) / (stddev(i) * stddev(j));
     77        corr(j, i) = corr(i, j);
     78      }
     79    }
     80    return corr;
    4681  }
    4782
  • trunk/yat/statistics/utility.h

    r3550 r3745  
    3636#include "yat/utility/deprecate.h"
    3737#include "yat/utility/iterator_traits.h"
     38#include "yat/utility/Matrix.h"
    3839#include "yat/utility/sort_index.h"
    3940#include "yat/utility/Vector.h"
     
    136137  double cdf_hypergeometric_P(unsigned int k, unsigned int n1,
    137138                              unsigned int n2, unsigned int t);
     139
     140  /**
     141     Calculates the Pearson correlation between each pair of columns.
     142
     143     \return a symmetric matrix with correlations
     144
     145     \since New in yat 0.16
     146   */
     147  utility::Matrix correlation(utility::Matrix x);
    138148
    139149
Note: See TracChangeset for help on using the changeset viewer.