src/math/covariance.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 "covariance.h"
  20 #include <gl/xalloc.h>
  21 #include "moments.h"
  22 #include <gsl/gsl_matrix.h>
  23 #include <data/case.h>
  24 #include <data/variable.h>
  25 #include <libpspp/misc.h>
  26
  27 #define n_MOMENTS (MOMENT_VARIANCE + 1)
  28
  29
  30 struct covariance
  31 {
  32   /* The variables for which the covariance matrix is to be calculated */
  33   size_t n_vars;
  34   const struct variable **vars;
  35
  36   /* The weight variable (or NULL if none) */
  37   const struct variable *wv;
  38
  39   /* A set of matrices containing the 0th, 1st and 2nd moments */
  40   gsl_matrix **moments;
  41
  42   /* The class of missing values to exclude */
  43   enum mv_class exclude;
  44
  45   /* An array of doubles representing the covariance matrix.
  46      Only the top triangle is included, and no diagonals */
  47   double *cm;
  48   int n_cm;
  49 };
  50
  51
  52
  53 /* Return a matrix containing the M th moments.
  54    The matrix is of size  NxN where N is the number of variables.
  55    Each row represents the moments of a variable.
  56    In the absence of missing values, the columns of this matrix will
  57    be identical.  If missing values are involved, then element (i,j)
  58    is the moment of the i th variable, when paired with the j th variable.
  59  */
  60 const gsl_matrix *
  61 covariance_moments (const struct covariance *cov, int m)
  62 {
  63   return cov->moments[m];
  64 }
  65
  66
  67
  68 /* Create a covariance struct */
  69 struct covariance *
  70 covariance_create (size_t n_vars, const struct variable **vars,
  71                    const struct variable *weight, enum mv_class exclude)
  72 {
  73   size_t i;
  74   struct covariance *cov = xmalloc (sizeof *cov);
  75   cov->vars = xmalloc (sizeof *cov->vars * n_vars);
  76
  77   cov->wv = weight;
  78   cov->n_vars = n_vars;
  79
  80   for (i = 0; i < n_vars; ++i)
  81     cov->vars[i] = vars[i];
  82
  83   cov->moments = xmalloc (sizeof *cov->moments * n_MOMENTS);
  84
  85   for (i = 0; i < n_MOMENTS; ++i)
  86     cov->moments[i] = gsl_matrix_calloc (n_vars, n_vars);
  87
  88   cov->exclude = exclude;
  89
  90   cov->n_cm = (n_vars * (n_vars - 1)  ) / 2;
  91
  92   cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm);
  93
  94   return cov;
  95 }
  96
  97 /* Return an integer, which can be used to index
  98    into COV->cm, to obtain the I, J th element
  99    of the covariance matrix.  If COV->cm does not
 100    contain that element, then a negative value
 101    will be returned.
 102 */
 103 static int
 104 cm_idx (const struct covariance *cov, int i, int j)
 105 {
 106   int as;
 107   const int n2j = cov->n_vars - 2 - j;
 108   const int nj = cov->n_vars - 2 ;
 109
 110   assert (i >= 0);
 111   assert (j < cov->n_vars);
 112
 113   if ( i == 0)
 114     return -1;
 115
 116   if (j >= cov->n_vars - 1)
 117     return -1;
 118
 119   if ( i <= j)
 120     return -1 ;
 121
 122   as = nj * (nj + 1) ;
 123   as -= n2j * (n2j + 1) ;
 124   as /= 2;
 125
 126   return i - 1 + as;
 127 }
 128
 129
 130 /* Call this function for every case in the data set */
 131 void
 132 covariance_accumulate (struct covariance *cov, const struct ccase *c)
 133 {
 134   size_t i, j, m;
 135   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
 136
 137   for (i = 0 ; i < cov->n_vars; ++i)
 138     {
 139       const union value *val1 = case_data (c, cov->vars[i]);
 140
 141       if ( var_is_value_missing (cov->vars[i], val1, cov->exclude))
 142         continue;
 143
 144       for (j = 0 ; j < cov->n_vars; ++j)
 145         {
 146           double pwr = 1.0;
 147           int idx;
 148           const union value *val2 = case_data (c, cov->vars[j]);
 149
 150           if ( var_is_value_missing (cov->vars[j], val2, cov->exclude))
 151             continue;
 152
 153           idx = cm_idx (cov, i, j);
 154           if (idx >= 0)
 155             {
 156               cov->cm [idx] += val1->f * val2->f;
 157             }
 158
 159           for (m = 0 ; m < n_MOMENTS; ++m)
 160             {
 161               double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 162
 163               *x += pwr * weight;
 164               pwr *= val1->f;
 165             }
 166         }
 167     }
 168 }
 169
 170
 171 /*
 172    Allocate and return a gsl_matrix containing the covariances of the
 173    data.
 174 */
 175 static gsl_matrix *
 176 cm_to_gsl (struct covariance *cov)
 177 {
 178   int i, j;
 179   gsl_matrix *m = gsl_matrix_calloc (cov->n_vars, cov->n_vars);
 180
 181   /* Copy the non-diagonal elements from cov->cm */
 182   for ( j = 0 ; j < cov->n_vars - 1; ++j)
 183     {
 184       for (i = j+1 ; i < cov->n_vars; ++i)
 185         {
 186           double x = cov->cm [cm_idx (cov, i, j)];
 187           gsl_matrix_set (m, i, j, x);
 188           gsl_matrix_set (m, j, i, x);
 189         }
 190     }
 191
 192   /* Copy the diagonal elements from cov->moments[2] */
 193   for (j = 0 ; j < cov->n_vars ; ++j)
 194     {
 195       double sigma = gsl_matrix_get (cov->moments[2], j, j);
 196       gsl_matrix_set (m, j, j, sigma);
 197     }
 198
 199   return m;
 200 }
 201
 202
 203
 204 /*
 205    Return a pointer to gsl_matrix containing the pairwise covariances.
 206    The matrix remains owned by the COV object, and must not be freed.
 207    Call this function only after all data have been accumulated.
 208 */
 209 const gsl_matrix *
 210 covariance_calculate (struct covariance *cov)
 211 {
 212   size_t i, j;
 213   size_t m;
 214
 215   for (m = 0; m < n_MOMENTS; ++m)
 216     {
 217       /* Divide the moments by the number of samples */
 218       if ( m > 0)
 219         {
 220           for (i = 0 ; i < cov->n_vars; ++i)
 221             {
 222               for (j = 0 ; j < cov->n_vars; ++j)
 223                 {
 224                   double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 225                   *x /= gsl_matrix_get (cov->moments[0], i, j);
 226
 227                   if ( m == MOMENT_VARIANCE)
 228                     *x -= pow2 (gsl_matrix_get (cov->moments[1], i, j));
 229                 }
 230             }
 231         }
 232     }
 233
 234   /* Centre the moments */
 235   for ( j = 0 ; j < cov->n_vars - 1; ++j)
 236     {
 237       for (i = j + 1 ; i < cov->n_vars; ++i)
 238         {
 239           double *x = &cov->cm [cm_idx (cov, i, j)];
 240
 241           *x /= gsl_matrix_get (cov->moments[0], i, j);
 242
 243           *x -=
 244             gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)
 245             *
 246             gsl_matrix_get (cov->moments[MOMENT_MEAN], j, i);
 247         }
 248     }
 249
 250   return cm_to_gsl (cov);
 251 }
 252
 253
 254 /* Destroy the COV object */
 255 void
 256 covariance_destroy (struct covariance *cov)
 257 {
 258   size_t i;
 259   free (cov->vars);
 260
 261   for (i = 0; i < n_MOMENTS; ++i)
 262     gsl_matrix_free (cov->moments[i]);
 263
 264   free (cov->moments);
 265   free (cov->cm);
 266   free (cov);
 267 }