1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2009, 2011 Free Software Foundation, Inc.
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.
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.
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/>. */
20 #include "math/correlation.h"
22 #include <gsl/gsl_matrix.h>
23 #include <gsl/gsl_cdf.h>
26 #include "libpspp/misc.h"
28 #include "gl/minmax.h"
32 significance_of_correlation (double rho, double w)
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));
45 return gsl_cdf_tdist_Q (t, w - 2);
47 return gsl_cdf_tdist_P (t, w - 2);
51 correlation_from_covariance (const gsl_matrix *cv, const gsl_matrix *v)
54 gsl_matrix *corr = gsl_matrix_calloc (cv->size1, cv->size2);
56 for (i = 0 ; i < cv->size1; ++i)
58 for (j = 0 ; j < cv->size2; ++j)
60 double rho = gsl_matrix_get (cv, i, j);
62 rho /= sqrt (gsl_matrix_get (v, i, j))
64 sqrt (gsl_matrix_get (v, j, i));
66 gsl_matrix_set (corr, i, j, rho);
74 covariance_from_correlation (const gsl_matrix *corr, const gsl_matrix *v)
77 assert (corr->size1 == corr->size2);
79 gsl_matrix *output = gsl_matrix_calloc (corr->size1, corr->size2);
81 for (i = 0 ; i < corr->size1; ++i)
83 for (j = 0 ; j < corr->size2; ++j)
85 double r = gsl_matrix_get (corr, i, j);
87 r *= sqrt (gsl_matrix_get (v, i, j))
89 sqrt (gsl_matrix_get (v, j, i));
91 gsl_matrix_set (output, i, j, r);