src/math/correlation.c

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