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

add function calculating correlation matrix

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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
Note: See TracChangeset for help on using the changeset viewer.