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 <libpspp/assertion.h>
  20 #include "covariance.h"
  21 #include <gl/xalloc.h>
  22 #include "moments.h"
  23 #include <gsl/gsl_matrix.h>
  24 #include <data/case.h>
  25 #include <data/variable.h>
  26 #include <libpspp/misc.h>
  27 #include "categoricals.h"
  28
  29 #define n_MOMENTS (MOMENT_VARIANCE + 1)
  30
  31
  32 struct covariance
  33 {
  34   /* The variables for which the covariance matrix is to be calculated. */
  35   size_t n_vars;
  36   const struct variable **vars;
  37
  38   /* Categorical variables. */
  39   struct categoricals *categoricals;
  40
  41   /* Array containing number of categories per categorical variable. */
  42   size_t *n_categories;
  43
  44   /* Dimension of the covariance matrix. */
  45   size_t dim;
  46
  47   /* The weight variable (or NULL if none) */
  48   const struct variable *wv;
  49
  50   /* A set of matrices containing the 0th, 1st and 2nd moments */
  51   gsl_matrix **moments;
  52
  53   /* The class of missing values to exclude */
  54   enum mv_class exclude;
  55
  56   /* An array of doubles representing the covariance matrix.
  57      Only the top triangle is included, and no diagonals */
  58   double *cm;
  59   int n_cm;
  60
  61   /* 1 for single pass algorithm;
  62      2 for double pass algorithm
  63   */
  64   short passes;
  65
  66   /*
  67     0 : No pass has  been made
  68     1 : First pass has been started
  69     2 : Second pass has been
  70
  71     IE: How many passes have been (partially) made. */
  72   short state;
  73
  74   /* Flags indicating that the first case has been seen */
  75   bool pass_one_first_case_seen;
  76   bool pass_two_first_case_seen;
  77 };
  78
  79
  80
  81 /* Return a matrix containing the M th moments.
  82    The matrix is of size  NxN where N is the number of variables.
  83    Each row represents the moments of a variable.
  84    In the absence of missing values, the columns of this matrix will
  85    be identical.  If missing values are involved, then element (i,j)
  86    is the moment of the i th variable, when paired with the j th variable.
  87  */
  88 const gsl_matrix *
  89 covariance_moments (const struct covariance *cov, int m)
  90 {
  91   return cov->moments[m];
  92 }
  93
  94
  95
  96 /* Create a covariance struct.
  97  */
  98 struct covariance *
  99 covariance_1pass_create (size_t n_vars, const struct variable **vars,
 100                          const struct variable *weight, enum mv_class exclude)
 101 {
 102   size_t i;
 103   struct covariance *cov = xmalloc (sizeof *cov);
 104
 105   cov->passes = 1;
 106   cov->state = 0;
 107   cov->pass_one_first_case_seen = cov->pass_two_first_case_seen = false;
 108
 109   cov->vars = vars;
 110
 111   cov->wv = weight;
 112   cov->n_vars = n_vars;
 113   cov->dim = n_vars;
 114
 115   cov->moments = xmalloc (sizeof *cov->moments * n_MOMENTS);
 116
 117   for (i = 0; i < n_MOMENTS; ++i)
 118     cov->moments[i] = gsl_matrix_calloc (n_vars, n_vars);
 119
 120   cov->exclude = exclude;
 121
 122   cov->n_cm = (n_vars * (n_vars - 1)  ) / 2;
 123
 124   cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm);
 125
 126   return cov;
 127 }
 128
 129 /*
 130   Create a covariance struct for a two-pass algorithm. If categorical
 131   variables are involed, the dimension cannot be know until after the
 132   first data pass, so the actual covariances will not be allocated
 133   until then.
 134  */
 135 struct covariance *
 136 covariance_2pass_create (size_t n_vars, const struct variable **vars,
 137                          size_t n_catvars, const struct variable **catvars,
 138                          const struct variable *wv, enum mv_class exclude)
 139 {
 140   size_t i;
 141   struct covariance *cov = xmalloc (sizeof *cov);
 142
 143   cov->passes = 2;
 144   cov->state = 0;
 145   cov->pass_one_first_case_seen = cov->pass_two_first_case_seen = false;
 146
 147   cov->vars = vars;
 148
 149   cov->wv = wv;
 150   cov->n_vars = n_vars;
 151   cov->dim = n_vars;
 152
 153   cov->moments = xmalloc (sizeof *cov->moments * n_MOMENTS);
 154
 155   for (i = 0; i < n_MOMENTS; ++i)
 156     cov->moments[i] = gsl_matrix_calloc (n_vars, n_vars);
 157
 158   cov->exclude = exclude;
 159
 160   cov->n_cm = - 1;
 161   cov->cm = NULL;
 162
 163   cov->categoricals = categoricals_create (catvars, n_catvars, wv);
 164
 165   return cov;
 166 }
 167
 168 /* Return an integer, which can be used to index
 169    into COV->cm, to obtain the I, J th element
 170    of the covariance matrix.  If COV->cm does not
 171    contain that element, then a negative value
 172    will be returned.
 173 */
 174 static int
 175 cm_idx (const struct covariance *cov, int i, int j)
 176 {
 177   int as;
 178   const int n2j = cov->n_vars - 2 - j;
 179   const int nj = cov->n_vars - 2 ;
 180
 181   assert (i >= 0);
 182   assert (j < cov->n_vars);
 183
 184   if ( i == 0)
 185     return -1;
 186
 187   if (j >= cov->n_vars - 1)
 188     return -1;
 189
 190   if ( i <= j)
 191     return -1 ;
 192
 193   as = nj * (nj + 1) ;
 194   as -= n2j * (n2j + 1) ;
 195   as /= 2;
 196
 197   return i - 1 + as;
 198 }
 199
 200 static void
 201 dump_matrix (const gsl_matrix *m)
 202 {
 203   size_t i, j;
 204
 205   for (i = 0 ; i < m->size1; ++i)
 206     {
 207       for (j = 0 ; j < m->size2; ++j)
 208         printf ("%02f ", gsl_matrix_get (m, i, j));
 209       printf ("\n");
 210     }
 211 }
 212
 213 /* Call this function for every case in the data set */
 214 void
 215 covariance_accumulate_pass1 (struct covariance *cov, const struct ccase *c)
 216 {
 217   size_t i, j, m;
 218   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
 219
 220   assert (cov->passes == 2);
 221   if (!cov->pass_one_first_case_seen)
 222     {
 223       assert (cov->state == 0);
 224       cov->state = 1;
 225     }
 226
 227   categoricals_update (cov->categoricals, c);
 228
 229   for (i = 0 ; i < cov->n_vars; ++i)
 230     {
 231       const union value *val1 = case_data (c, cov->vars[i]);
 232
 233       if ( var_is_value_missing (cov->vars[i], val1, cov->exclude))
 234         continue;
 235
 236       for (j = 0 ; j < cov->n_vars; ++j)
 237         {
 238           double pwr = 1.0;
 239           const union value *val2 = case_data (c, cov->vars[j]);
 240
 241           if ( var_is_value_missing (cov->vars[j], val2, cov->exclude))
 242             continue;
 243
 244           for (m = 0 ; m <= MOMENT_MEAN; ++m)
 245             {
 246               double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 247
 248               *x += pwr * weight;
 249               pwr *= val1->f;
 250             }
 251         }
 252     }
 253
 254   cov->pass_one_first_case_seen = true;
 255 }
 256
 257
 258 /* Call this function for every case in the data set */
 259 void
 260 covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
 261 {
 262   size_t i, j;
 263   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
 264
 265   assert (cov->passes == 2);
 266   assert (cov->state >= 1);
 267
 268   if (! cov->pass_two_first_case_seen)
 269     {
 270       assert (cov->state == 1);
 271       cov->state = 2;
 272
 273       cov->dim = cov->n_vars + categoricals_total (cov->categoricals);
 274       cov->n_cm = (cov->dim * (cov->dim - 1)  ) / 2;
 275       cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm);
 276
 277       /* Divide the means by the number of samples */
 278       for (i = 0; i < cov->n_vars; ++i)
 279         {
 280           for (j = 0; j < cov->n_vars; ++j)
 281             {
 282               double *x = gsl_matrix_ptr (cov->moments[MOMENT_MEAN], i, j);
 283               *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 284             }
 285         }
 286     }
 287
 288   for (i = 0 ; i < cov->n_vars; ++i)
 289     {
 290       const union value *val1 = case_data (c, cov->vars[i]);
 291
 292       if ( var_is_value_missing (cov->vars[i], val1, cov->exclude))
 293         continue;
 294
 295       for (j = 0 ; j < cov->n_vars; ++j)
 296         {
 297           int idx;
 298           double ss ;
 299           const union value *val2 = case_data (c, cov->vars[j]);
 300
 301           const double s = pow2 (val1->f - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)) * weight;
 302
 303           if ( var_is_value_missing (cov->vars[j], val2, cov->exclude))
 304             continue;
 305
 306           {
 307             double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
 308             *x += s;
 309           }
 310
 311           ss =
 312             (val1->f - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
 313             *
 314             (val2->f - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
 315             * weight
 316             ;
 317
 318           idx = cm_idx (cov, i, j);
 319           if (idx >= 0)
 320             {
 321               cov->cm [idx] += ss;
 322             }
 323
 324         }
 325     }
 326
 327   cov->pass_two_first_case_seen = true;
 328 }
 329
 330
 331 /* Call this function for every case in the data set.
 332    After all cases have been passed, call covariance_calculate
 333  */
 334 void
 335 covariance_accumulate (struct covariance *cov, const struct ccase *c)
 336 {
 337   size_t i, j, m;
 338   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
 339
 340   assert (cov->passes == 1);
 341
 342   if ( !cov->pass_one_first_case_seen)
 343     {
 344       assert ( cov->state == 0);
 345       cov->state = 1;
 346     }
 347
 348   for (i = 0 ; i < cov->n_vars; ++i)
 349     {
 350       const union value *val1 = case_data (c, cov->vars[i]);
 351
 352       if ( var_is_value_missing (cov->vars[i], val1, cov->exclude))
 353         continue;
 354
 355       for (j = 0 ; j < cov->n_vars; ++j)
 356         {
 357           double pwr = 1.0;
 358           int idx;
 359           const union value *val2 = case_data (c, cov->vars[j]);
 360
 361           if ( var_is_value_missing (cov->vars[j], val2, cov->exclude))
 362             continue;
 363
 364           idx = cm_idx (cov, i, j);
 365           if (idx >= 0)
 366             {
 367               cov->cm [idx] += val1->f * val2->f * weight;
 368             }
 369
 370           for (m = 0 ; m < n_MOMENTS; ++m)
 371             {
 372               double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 373
 374               *x += pwr * weight;
 375               pwr *= val1->f;
 376             }
 377         }
 378     }
 379
 380   cov->pass_one_first_case_seen = true;
 381 }
 382
 383
 384 /*
 385    Allocate and return a gsl_matrix containing the covariances of the
 386    data.
 387 */
 388 static gsl_matrix *
 389 cm_to_gsl (struct covariance *cov)
 390 {
 391   int i, j;
 392   gsl_matrix *m = gsl_matrix_calloc (cov->n_vars, cov->n_vars);
 393
 394   /* Copy the non-diagonal elements from cov->cm */
 395   for ( j = 0 ; j < cov->n_vars - 1; ++j)
 396     {
 397       for (i = j+1 ; i < cov->n_vars; ++i)
 398         {
 399           double x = cov->cm [cm_idx (cov, i, j)];
 400           gsl_matrix_set (m, i, j, x);
 401           gsl_matrix_set (m, j, i, x);
 402         }
 403     }
 404
 405   /* Copy the diagonal elements from cov->moments[2] */
 406   for (j = 0 ; j < cov->n_vars ; ++j)
 407     {
 408       double sigma = gsl_matrix_get (cov->moments[2], j, j);
 409       gsl_matrix_set (m, j, j, sigma);
 410     }
 411
 412   return m;
 413 }
 414
 415
 416 static const gsl_matrix *
 417 covariance_calculate_double_pass (struct covariance *cov)
 418 {
 419   size_t i, j;
 420   for (i = 0 ; i < cov->n_vars; ++i)
 421     {
 422       for (j = 0 ; j < cov->n_vars; ++j)
 423         {
 424           int idx;
 425           double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
 426           *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 427
 428           idx = cm_idx (cov, i, j);
 429           if ( idx >= 0)
 430             {
 431               x = &cov->cm [idx];
 432               *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 433             }
 434         }
 435     }
 436
 437   return  cm_to_gsl (cov);
 438 }
 439
 440 static const gsl_matrix *
 441 covariance_calculate_single_pass (struct covariance *cov)
 442 {
 443   size_t i, j;
 444   size_t m;
 445
 446   for (m = 0; m < n_MOMENTS; ++m)
 447     {
 448       /* Divide the moments by the number of samples */
 449       if ( m > 0)
 450         {
 451           for (i = 0 ; i < cov->n_vars; ++i)
 452             {
 453               for (j = 0 ; j < cov->n_vars; ++j)
 454                 {
 455                   double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 456                   *x /= gsl_matrix_get (cov->moments[0], i, j);
 457
 458                   if ( m == MOMENT_VARIANCE)
 459                     *x -= pow2 (gsl_matrix_get (cov->moments[1], i, j));
 460                 }
 461             }
 462         }
 463     }
 464
 465   /* Centre the moments */
 466   for ( j = 0 ; j < cov->n_vars - 1; ++j)
 467     {
 468       for (i = j + 1 ; i < cov->n_vars; ++i)
 469         {
 470           double *x = &cov->cm [cm_idx (cov, i, j)];
 471
 472           *x /= gsl_matrix_get (cov->moments[0], i, j);
 473
 474           *x -=
 475             gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)
 476             *
 477             gsl_matrix_get (cov->moments[MOMENT_MEAN], j, i);
 478         }
 479     }
 480
 481   return cm_to_gsl (cov);
 482 }
 483
 484
 485
 486 /*
 487    Return a pointer to gsl_matrix containing the pairwise covariances.
 488    The matrix remains owned by the COV object, and must not be freed.
 489    Call this function only after all data have been accumulated.
 490 */
 491 const gsl_matrix *
 492 covariance_calculate (struct covariance *cov)
 493 {
 494   assert ( cov->state > 0 );
 495
 496   switch (cov->passes)
 497     {
 498     case 1:
 499       return covariance_calculate_single_pass (cov);
 500       break;
 501     case 2:
 502       return covariance_calculate_double_pass (cov);
 503       break;
 504     default:
 505       NOT_REACHED ();
 506     }
 507 }
 508
 509
 510
 511
 512 /* Destroy the COV object */
 513 void
 514 covariance_destroy (struct covariance *cov)
 515 {
 516   size_t i;
 517   free (cov->vars);
 518
 519   for (i = 0; i < n_MOMENTS; ++i)
 520     gsl_matrix_free (cov->moments[i]);
 521
 522   free (cov->moments);
 523   free (cov->cm);
 524   free (cov);
 525 }