pintos-os.org Git - pspp/blob - src/math/correlation.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2009, 2011 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 #include <config.h>
  18
  19 #include <assert.h>
  20 #include "math/correlation.h"
  21
  22 #include <gsl/gsl_matrix.h>
  23 #include <gsl/gsl_cdf.h>
  24 #include <math.h>
  25
  26 #include "libpspp/misc.h"
  27
  28 #include "gl/minmax.h"
  29
  30
  31 double
  32 significance_of_correlation (double rho, double w)
  33 {
  34   double t = w - 2;
  35
  36   /* |rho| will mathematically always be in the range [0, 1.0].  Inaccurate
  37      calculations sometimes cause it to be slightly greater than 1.0, so
  38      force it into the correct range to avoid NaN from sqrt(). */
  39   t /= 1 - MIN (1, pow2 (rho));
  40
  41   t = sqrt (t);
  42   t *= rho;
  43
  44   if (t > 0)
  45     return  gsl_cdf_tdist_Q (t, w - 2);
  46   else
  47     return  gsl_cdf_tdist_P (t, w - 2);
  48 }
  49
  50 gsl_matrix *
  51 correlation_from_covariance (const gsl_matrix *cv, const gsl_matrix *v)
  52 {
  53   size_t i, j;
  54   gsl_matrix *corr = gsl_matrix_calloc (cv->size1, cv->size2);
  55
  56   for (i = 0 ; i < cv->size1; ++i)
  57     {
  58       for (j = 0 ; j < cv->size2; ++j)
  59         {
  60           double rho = gsl_matrix_get (cv, i, j);
  61
  62           rho /= sqrt (gsl_matrix_get (v, i, j))
  63             *
  64             sqrt (gsl_matrix_get (v, j, i));
  65
  66           gsl_matrix_set (corr, i, j, rho);
  67         }
  68     }
  69
  70   return corr;
  71 }
  72
  73 gsl_matrix *
  74 covariance_from_correlation (const gsl_matrix *corr, const gsl_matrix *v)
  75 {
  76   size_t i, j;
  77   assert (corr->size1 == corr->size2);
  78
  79   gsl_matrix *output = gsl_matrix_calloc (corr->size1, corr->size2);
  80
  81   for (i = 0 ; i < corr->size1; ++i)
  82     {
  83       for (j = 0 ; j < corr->size2; ++j)
  84         {
  85           double r = gsl_matrix_get (corr, i, j);
  86
  87           r *= sqrt (gsl_matrix_get (v, i, j))
  88             *
  89             sqrt (gsl_matrix_get (v, j, i));
  90
  91           gsl_matrix_set (output, i, j, r);
  92         }
  93     }
  94
  95   return output;
  96 }