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