X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fmath%2Fcorrelation.c;h=bac81e4349360ef902a7e99b63768c0cafed42af;hb=b9f712f8911170fcc049c6d7d83b6887185de4dd;hp=303b48bcfaccb89098fabaa426f628d50323a2a3;hpb=50ac6802fc247814dc4dd6232f6304b928a2d78b;p=pspp

diff --git a/src/math/correlation.c b/src/math/correlation.c
index 303b48bcfa..bac81e4349 100644
--- a/src/math/correlation.c
+++ b/src/math/correlation.c
@@ -1,5 +1,5 @@
 /* PSPP - a program for statistical analysis.
-   Copyright (C) 2009 Free Software Foundation, Inc.
+   Copyright (C) 2009, 2011 Free Software Foundation, Inc.
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
@@ -16,23 +16,31 @@
 
 #include <config.h>
 
-#include "correlation.h"
+#include <assert.h>
+#include "math/correlation.h"
 
 #include <gsl/gsl_matrix.h>
 #include <gsl/gsl_cdf.h>
 #include <math.h>
-#include <libpspp/misc.h>
-#include "minmax.h"
+
+#include "libpspp/misc.h"
+
+#include "gl/minmax.h"
 
 
 double
 significance_of_correlation (double rho, double w)
 {
   double t = w - 2;
+
+  /* |rho| will mathematically always be in the range [0, 1.0].  Inaccurate
+     calculations sometimes cause it to be slightly greater than 1.0, so
+     force it into the correct range to avoid NaN from sqrt(). */
   t /= 1 - MIN (1, pow2 (rho));
+
   t = sqrt (t);
   t *= rho;
-  
+
   if (t > 0)
     return  gsl_cdf_tdist_Q (t, w - 2);
   else
@@ -44,20 +52,45 @@ correlation_from_covariance (const gsl_matrix *cv, const gsl_matrix *v)
 {
   size_t i, j;
   gsl_matrix *corr = gsl_matrix_calloc (cv->size1, cv->size2);
-  
+
   for (i = 0 ; i < cv->size1; ++i)
     {
       for (j = 0 ; j < cv->size2; ++j)
 	{
 	  double rho = gsl_matrix_get (cv, i, j);
-	  
+
 	  rho /= sqrt (gsl_matrix_get (v, i, j))
-	    * 
+	    *
 	    sqrt (gsl_matrix_get (v, j, i));
-	  
+
 	  gsl_matrix_set (corr, i, j, rho);
 	}
     }
-  
+
   return corr;
 }
+
+gsl_matrix *
+covariance_from_correlation (const gsl_matrix *corr, const gsl_matrix *v)
+{
+  size_t i, j;
+  assert (corr->size1 == corr->size2);
+
+  gsl_matrix *output = gsl_matrix_calloc (corr->size1, corr->size2);
+
+  for (i = 0 ; i < corr->size1; ++i)
+    {
+      for (j = 0 ; j < corr->size2; ++j)
+	{
+	  double r = gsl_matrix_get (corr, i, j);
+
+	  r *= sqrt (gsl_matrix_get (v, i, j))
+	    *
+	    sqrt (gsl_matrix_get (v, j, i));
+
+	  gsl_matrix_set (output, i, j, r);
+	}
+    }
+
+  return output;
+}