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 "math/correlation.h"
  20
  21 #include <gsl/gsl_matrix.h>
  22 #include <gsl/gsl_cdf.h>
  23 #include <math.h>
  24
  25 #include "libpspp/misc.h"
  26
  27 #include "gl/minmax.h"
  28
  29
  30 double
  31 significance_of_correlation (double rho, double w)
  32 {
  33   double t = w - 2;
  34
  35   /* |rho| will mathematically always be in the range [0, 1.0].  Inaccurate
  36      calculations sometimes cause it to be slightly greater than 1.0, so
  37      force it into the correct range to avoid NaN from sqrt(). */
  38   t /= 1 - MIN (1, pow2 (rho));
  39
  40   t = sqrt (t);
  41   t *= rho;
  42
  43   if (t > 0)
  44     return  gsl_cdf_tdist_Q (t, w - 2);
  45   else
  46     return  gsl_cdf_tdist_P (t, w - 2);
  47 }
  48
  49 gsl_matrix *
  50 correlation_from_covariance (const gsl_matrix *cv, const gsl_matrix *v)
  51 {
  52   size_t i, j;
  53   gsl_matrix *corr = gsl_matrix_calloc (cv->size1, cv->size2);
  54
  55   for (i = 0 ; i < cv->size1; ++i)
  56     {
  57       for (j = 0 ; j < cv->size2; ++j)
  58         {
  59           double rho = gsl_matrix_get (cv, i, j);
  60
  61           rho /= sqrt (gsl_matrix_get (v, i, j))
  62             *
  63             sqrt (gsl_matrix_get (v, j, i));
  64
  65           gsl_matrix_set (corr, i, j, rho);
  66         }
  67     }
  68
  69   return corr;
  70 }