src/math/covariance-matrix.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2008 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 /*
  18   Create and update the values in the covariance matrix.
  19 */
  20 #include <assert.h>
  21 #include <config.h>
  22 #include <data/variable.h>
  23 #include <data/value.h>
  24 #include "covariance-matrix.h"
  25 #include "moments.h"
  26
  27 /*
  28   The covariances are stored in a DESIGN_MATRIX structure.
  29  */
  30 struct design_matrix *
  31 covariance_matrix_create (int n_variables, const struct variable *v_variables[])
  32 {
  33   return design_matrix_create (n_variables, v_variables, (size_t) n_variables);
  34 }
  35
  36 void covariance_matrix_destroy (struct design_matrix *x)
  37 {
  38   design_matrix_destroy (x);
  39 }
  40
  41 /*
  42   Update the covariance matrix with the new entries, assuming that ROW
  43   corresponds to a categorical variable and V2 is numeric.
  44  */
  45 static void
  46 covariance_update_categorical_numeric (struct design_matrix *cov, double mean,
  47                           size_t row,
  48                           const struct variable *v2, double x, const union value *val2)
  49 {
  50   size_t col;
  51   double tmp;
  52
  53   assert (var_is_numeric (v2));
  54
  55   col = design_matrix_var_to_column (cov, v2);
  56   assert (val2 != NULL);
  57   tmp = gsl_matrix_get (cov->m, row, col);
  58   gsl_matrix_set (cov->m, row, col, (val2->f - mean) * x + tmp);
  59   gsl_matrix_set (cov->m, col, row, (val2->f - mean) * x + tmp);
  60 }
  61 static void
  62 column_iterate (struct design_matrix *cov, const struct variable *v,
  63                 double ssize, double x, const union value *val1, size_t row)
  64 {
  65   size_t col;
  66   size_t i;
  67   double y;
  68   double tmp;
  69   const union value *tmp_val;
  70
  71   col = design_matrix_var_to_column (cov, v);
  72   for (i = 0; i < cat_get_n_categories (v) - 1; i++)
  73     {
  74       col += i;
  75       y = -1.0 * cat_get_category_count (i, v) / ssize;
  76       tmp_val = cat_subscript_to_value (i, v);
  77       if (compare_values (tmp_val, val1, var_get_width (v)))
  78         {
  79           y += -1.0;
  80         }
  81       tmp = gsl_matrix_get (cov->m, row, col);
  82       gsl_matrix_set (cov->m, row, col, x * y + tmp);
  83       gsl_matrix_set (cov->m, col, row, x * y + tmp);
  84     }
  85 }
  86 /*
  87   Call this function in the second data pass. The central moments are
  88   MEAN1 and MEAN2. Any categorical variables should already have their
  89   values summarized in in its OBS_VALS element.
  90  */
  91 void covariance_pass_two (struct design_matrix *cov, double mean1, double mean2,
  92                           double ssize, const struct variable *v1,
  93                           const struct variable *v2, const union value *val1, const union value *val2)
  94 {
  95   size_t row;
  96   size_t col;
  97   size_t i;
  98   double x;
  99   const union value *tmp_val;
 100
 101   if (var_is_alpha (v1))
 102     {
 103       row = design_matrix_var_to_column (cov, v1);
 104       for (i = 0; i < cat_get_n_categories (v1) - 1; i++)
 105         {
 106           row += i;
 107           x = -1.0 * cat_get_category_count (i, v1) / ssize;
 108           tmp_val = cat_subscript_to_value (i, v1);
 109           if (compare_values (tmp_val, val1, var_get_width (v1)))
 110             {
 111               x += 1.0;
 112             }
 113           if (var_is_numeric (v2))
 114             {
 115               covariance_update_categorical_numeric (cov, mean2, row,
 116                                                      v2, x, val2);
 117             }
 118           else
 119             {
 120               column_iterate (cov, v1, ssize, x, val1, row);
 121               column_iterate (cov, v2, ssize, x, val2, row);
 122             }
 123         }
 124     }
 125   else if (var_is_alpha (v2))
 126     {
 127       /*
 128         Reverse the orders of V1, V2, etc. and put ourselves back
 129         in the previous IF scope.
 130        */
 131       covariance_pass_two (cov, mean2, mean1, ssize, v2, v1, val2, val1);
 132     }
 133   else
 134     {
 135       /*
 136         Both variables are numeric.
 137       */
 138       row = design_matrix_var_to_column (cov, v1);
 139       col = design_matrix_var_to_column (cov, v2);
 140       x = (val1->f - mean1) * (val2->f - mean2);
 141       x += gsl_matrix_get (cov->m, col, row);
 142       gsl_matrix_set (cov->m, row, col, x);
 143       gsl_matrix_set (cov->m, col, row, x);
 144     }
 145 }
 146