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