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 V1
  43   is categorical and V2 is numeric.
  44  */
  45 static void
  46 covariance_update_categorical_numeric (struct design_matrix *cov, double mean,
  47                           double weight, double ssize, const struct variable *v1,
  48                           const struct variable *v2, const union value *val1, const union value *val2)
  49 {
  50   double x;
  51   size_t i;
  52   size_t col;
  53   size_t row;
  54
  55   assert (var_is_alpha (v1));
  56   assert (var_is_numeric (v2));
  57
  58   row = design_matrix_var_to_column (cov, v1);
  59   col = design_matrix_var_to_column (cov, v2);
  60   for (i = 0; i < cat_get_n_categories (v1); i++)
  61     {
  62       row += i;
  63       x = -1.0 * cat_get_n_categories (v1) / ssize;
  64       if (i == cat_value_find (v1, val1))
  65         {
  66           x += 1.0;
  67         }
  68       assert (val2 != NULL);
  69       gsl_matrix_set (cov->m, row, col, (val2->f - mean) * x * weight);
  70     }
  71 }
  72
  73 /*
  74   Call this function in the first data pass. The central moments are
  75   MEAN1 and MEAN2. Any categorical variables should already have their
  76   values summarized in in its OBS_VALS element.
  77  */
  78 void covariance_pass_one (struct design_matrix *cov, double mean1, double mean2,
  79                           double weight, double ssize, const struct variable *v1,
  80                           const struct variable *v2, const union value *val1, const union value *val2)
  81 {
  82   size_t row;
  83   size_t col;
  84   size_t i;
  85   double x;
  86   double y;
  87
  88   if (var_is_alpha (v1))
  89     {
  90       if (var_is_numeric (v2))
  91         {
  92           covariance_update_categorical_numeric (cov, mean2, weight, ssize, v1,
  93                                                  v2, val1, val2);
  94         }
  95       else
  96         {
  97           row = design_matrix_var_to_column (cov, v1);
  98           col = design_matrix_var_to_column (cov, v2);
  99           for (i = 0; i < cat_get_n_categories (v2); i++)
 100             {
 101               col += i;
 102               y = -1.0 * cat_get_n_categories (v2) / ssize;
 103               if (i == cat_value_find (v2, val2))
 104                 {
 105                   y += 1.0;
 106                 }
 107               gsl_matrix_set (cov->m, row, col, x * y * weight);
 108               gsl_matrix_set (cov->m, col, row, x * y * weight);
 109             }
 110         }
 111     }
 112   else if (var_is_alpha (v2))
 113     {
 114       covariance_update_categorical_numeric (cov, mean1, weight, ssize, v2,
 115                                              v1, val2, val1);
 116     }
 117   else
 118     {
 119       /*
 120         Both variables are numeric.
 121       */
 122       row = design_matrix_var_to_column (cov, v1);
 123       col = design_matrix_var_to_column (cov, v2);
 124       x = (val1->f - mean1) * (val2->f - mean2) * weight;
 125       gsl_matrix_set (cov->m, row, col, x);
 126       gsl_matrix_set (cov->m, col, row, x);
 127     }
 128 }
 129