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 = xzalloc (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   if (cov->n_cm > 0)
 160     cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm);
 161   cov->categoricals = NULL;
 162
 163   return cov;
 164 }
 165
 166 /*
 167   Create a covariance struct for a two-pass algorithm. If categorical
 168   variables are involed, the dimension cannot be know until after the
 169   first data pass, so the actual covariances will not be allocated
 170   until then.
 171  */
 172 struct covariance *
 173 covariance_2pass_create (size_t n_vars, const struct variable **vars,
 174                          size_t n_catvars, const struct variable **catvars,
 175                          const struct variable *wv, enum mv_class exclude)
 176 {
 177   size_t i;
 178   struct covariance *cov = xmalloc (sizeof *cov);
 179
 180   cov->passes = 2;
 181   cov->state = 0;
 182   cov->pass_one_first_case_seen = cov->pass_two_first_case_seen = false;
 183
 184   cov->vars = vars;
 185
 186   cov->wv = wv;
 187   cov->n_vars = n_vars;
 188   cov->dim = n_vars;
 189
 190   cov->moments = xmalloc (sizeof *cov->moments * n_MOMENTS);
 191
 192   for (i = 0; i < n_MOMENTS; ++i)
 193     cov->moments[i] = gsl_matrix_calloc (n_vars, n_vars);
 194
 195   cov->exclude = exclude;
 196
 197   cov->n_cm = -1;
 198   cov->cm = NULL;
 199
 200   cov->categoricals = categoricals_create (catvars, n_catvars, wv, exclude);
 201
 202   return cov;
 203 }
 204
 205 /* Return an integer, which can be used to index
 206    into COV->cm, to obtain the I, J th element
 207    of the covariance matrix.  If COV->cm does not
 208    contain that element, then a negative value
 209    will be returned.
 210 */
 211 static int
 212 cm_idx (const struct covariance *cov, int i, int j)
 213 {
 214   int as;
 215   const int n2j = cov->dim - 2 - j;
 216   const int nj = cov->dim - 2 ;
 217
 218   assert (i >= 0);
 219   assert (j < cov->dim);
 220
 221   if ( i == 0)
 222     return -1;
 223
 224   if (j >= cov->dim - 1)
 225     return -1;
 226
 227   if ( i <= j)
 228     return -1 ;
 229
 230   as = nj * (nj + 1) ;
 231   as -= n2j * (n2j + 1) ;
 232   as /= 2;
 233
 234   return i - 1 + as;
 235 }
 236
 237
 238 /*
 239   Returns true iff the variable corresponding to the Ith element of the covariance matrix
 240    has a missing value for case C
 241 */
 242 static bool
 243 is_missing (const struct covariance *cov, int i, const struct ccase *c)
 244 {
 245   const struct variable *var = i < cov->n_vars ?
 246     cov->vars[i] :
 247     categoricals_get_variable_by_subscript (cov->categoricals, i - cov->n_vars);
 248
 249   const union value *val = case_data (c, var);
 250
 251   return var_is_value_missing (var, val, cov->exclude);
 252 }
 253
 254
 255 static double
 256 get_val (const struct covariance *cov, int i, const struct ccase *c)
 257 {
 258   if ( i < cov->n_vars)
 259     {
 260       const struct variable *var = cov->vars[i];
 261
 262       const union value *val = case_data (c, var);
 263
 264       return val->f;
 265     }
 266
 267   return categoricals_get_binary_by_subscript (cov->categoricals, i - cov->n_vars, c);
 268 }
 269
 270 void
 271 dump_matrix (const gsl_matrix *m)
 272 {
 273   size_t i, j;
 274
 275   for (i = 0 ; i < m->size1; ++i)
 276     {
 277       for (j = 0 ; j < m->size2; ++j)
 278         printf ("%02f ", gsl_matrix_get (m, i, j));
 279       printf ("\n");
 280     }
 281 }
 282
 283 /* Call this function for every case in the data set */
 284 void
 285 covariance_accumulate_pass1 (struct covariance *cov, const struct ccase *c)
 286 {
 287   size_t i, j, m;
 288   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
 289
 290   assert (cov->passes == 2);
 291   if (!cov->pass_one_first_case_seen)
 292     {
 293       assert (cov->state == 0);
 294       cov->state = 1;
 295     }
 296
 297   categoricals_update (cov->categoricals, c);
 298
 299   for (i = 0 ; i < cov->dim; ++i)
 300     {
 301       double v1 = get_val (cov, i, c);
 302
 303       if ( is_missing (cov, i, c))
 304         continue;
 305
 306       for (j = 0 ; j < cov->dim; ++j)
 307         {
 308           double pwr = 1.0;
 309
 310           if ( is_missing (cov, j, c))
 311             continue;
 312
 313           for (m = 0 ; m <= MOMENT_MEAN; ++m)
 314             {
 315               double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 316
 317               *x += pwr * weight;
 318               pwr *= v1;
 319             }
 320         }
 321     }
 322
 323   cov->pass_one_first_case_seen = true;
 324 }
 325
 326
 327 /* Call this function for every case in the data set */
 328 void
 329 covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
 330 {
 331   size_t i, j;
 332   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
 333
 334   assert (cov->passes == 2);
 335   assert (cov->state >= 1);
 336
 337   if (! cov->pass_two_first_case_seen)
 338     {
 339       size_t m;
 340       assert (cov->state == 1);
 341       cov->state = 2;
 342
 343       cov->dim = cov->n_vars +
 344         categoricals_total (cov->categoricals) - categoricals_get_n_variables (cov->categoricals);
 345
 346       cov->n_cm = (cov->dim * (cov->dim - 1)  ) / 2;
 347       cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm);
 348
 349       /* Grow the moment matrices so that they're large enough to accommodate the
 350          categorical elements */
 351       for (i = 0; i < n_MOMENTS; ++i)
 352         {
 353           cov->moments[i] = resize_matrix (cov->moments[i], cov->dim);
 354         }
 355
 356       categoricals_done (cov->categoricals);
 357
 358       /* Populate the moments matrices with the categorical value elements */
 359       for (i = cov->n_vars; i < cov->dim; ++i)
 360         {
 361           for (j = 0 ; j < cov->dim ; ++j) /* FIXME: This is WRONG !!! */
 362             {
 363               double w = categoricals_get_weight_by_subscript (cov->categoricals, i - cov->n_vars);
 364
 365               gsl_matrix_set (cov->moments[MOMENT_NONE], i, j, w);
 366
 367               w = categoricals_get_sum_by_subscript (cov->categoricals, i - cov->n_vars);
 368
 369               gsl_matrix_set (cov->moments[MOMENT_MEAN], i, j, w);
 370             }
 371         }
 372
 373       /* FIXME: This is WRONG!!  It must be fixed to properly handle missing values.  For
 374        now it assumes there are none */
 375       for (m = 0 ; m < n_MOMENTS; ++m)
 376         {
 377           for (i = 0 ; i < cov->dim ; ++i)
 378             {
 379               double x = gsl_matrix_get (cov->moments[m], i, cov->n_vars -1);
 380               for (j = cov->n_vars; j < cov->dim; ++j)
 381                 {
 382                   gsl_matrix_set (cov->moments[m], i, j, x);
 383                 }
 384             }
 385         }
 386
 387       /* Divide the means by the number of samples */
 388       for (i = 0; i < cov->dim; ++i)
 389         {
 390           for (j = 0; j < cov->dim; ++j)
 391             {
 392               double *x = gsl_matrix_ptr (cov->moments[MOMENT_MEAN], i, j);
 393               *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 394             }
 395         }
 396     }
 397
 398   for (i = 0 ; i < cov->dim; ++i)
 399     {
 400       double v1 = get_val (cov, i, c);
 401
 402       if ( is_missing (cov, i, c))
 403         continue;
 404
 405       for (j = 0 ; j < cov->dim; ++j)
 406         {
 407           int idx;
 408           double ss ;
 409           double v2 = get_val (cov, j, c);
 410
 411           const double s = pow2 (v1 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)) * weight;
 412
 413           if ( is_missing (cov, j, c))
 414             continue;
 415
 416           {
 417             double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
 418             *x += s;
 419           }
 420
 421           ss =
 422             (v1 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
 423             *
 424             (v2 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
 425             * weight
 426             ;
 427
 428           idx = cm_idx (cov, i, j);
 429           if (idx >= 0)
 430             {
 431               cov->cm [idx] += ss;
 432             }
 433         }
 434     }
 435
 436   cov->pass_two_first_case_seen = true;
 437 }
 438
 439
 440 /* Call this function for every case in the data set.
 441    After all cases have been passed, call covariance_calculate
 442  */
 443 void
 444 covariance_accumulate (struct covariance *cov, const struct ccase *c)
 445 {
 446   size_t i, j, m;
 447   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
 448
 449   assert (cov->passes == 1);
 450
 451   if ( !cov->pass_one_first_case_seen)
 452     {
 453       assert ( cov->state == 0);
 454       cov->state = 1;
 455     }
 456
 457   for (i = 0 ; i < cov->dim; ++i)
 458     {
 459       const union value *val1 = case_data (c, cov->vars[i]);
 460
 461       if ( is_missing (cov, i, c))
 462         continue;
 463
 464       for (j = 0 ; j < cov->dim; ++j)
 465         {
 466           double pwr = 1.0;
 467           int idx;
 468           const union value *val2 = case_data (c, cov->vars[j]);
 469
 470           if ( is_missing (cov, j, c))
 471             continue;
 472
 473           idx = cm_idx (cov, i, j);
 474           if (idx >= 0)
 475             {
 476               cov->cm [idx] += val1->f * val2->f * weight;
 477             }
 478
 479           for (m = 0 ; m < n_MOMENTS; ++m)
 480             {
 481               double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 482
 483               *x += pwr * weight;
 484               pwr *= val1->f;
 485             }
 486         }
 487     }
 488
 489   cov->pass_one_first_case_seen = true;
 490 }
 491
 492
 493 /*
 494    Allocate and return a gsl_matrix containing the covariances of the
 495    data.
 496 */
 497 static gsl_matrix *
 498 cm_to_gsl (struct covariance *cov)
 499 {
 500   int i, j;
 501   gsl_matrix *m = gsl_matrix_calloc (cov->dim, cov->dim);
 502
 503   /* Copy the non-diagonal elements from cov->cm */
 504   for ( j = 0 ; j < cov->dim - 1; ++j)
 505     {
 506       for (i = j+1 ; i < cov->dim; ++i)
 507         {
 508           double x = cov->cm [cm_idx (cov, i, j)];
 509           gsl_matrix_set (m, i, j, x);
 510           gsl_matrix_set (m, j, i, x);
 511         }
 512     }
 513
 514   /* Copy the diagonal elements from cov->moments[2] */
 515   for (j = 0 ; j < cov->dim ; ++j)
 516     {
 517       double sigma = gsl_matrix_get (cov->moments[2], j, j);
 518       gsl_matrix_set (m, j, j, sigma);
 519     }
 520
 521   return m;
 522 }
 523
 524
 525 static const gsl_matrix *
 526 covariance_calculate_double_pass (struct covariance *cov)
 527 {
 528   size_t i, j;
 529   for (i = 0 ; i < cov->dim; ++i)
 530     {
 531       for (j = 0 ; j < cov->dim; ++j)
 532         {
 533           int idx;
 534           double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
 535           *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 536
 537           idx = cm_idx (cov, i, j);
 538           if ( idx >= 0)
 539             {
 540               x = &cov->cm [idx];
 541               *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 542             }
 543         }
 544     }
 545
 546   return  cm_to_gsl (cov);
 547 }
 548
 549 static const gsl_matrix *
 550 covariance_calculate_single_pass (struct covariance *cov)
 551 {
 552   size_t i, j;
 553   size_t m;
 554
 555   for (m = 0; m < n_MOMENTS; ++m)
 556     {
 557       /* Divide the moments by the number of samples */
 558       if ( m > 0)
 559         {
 560           for (i = 0 ; i < cov->dim; ++i)
 561             {
 562               for (j = 0 ; j < cov->dim; ++j)
 563                 {
 564                   double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 565                   *x /= gsl_matrix_get (cov->moments[0], i, j);
 566
 567                   if ( m == MOMENT_VARIANCE)
 568                     *x -= pow2 (gsl_matrix_get (cov->moments[1], i, j));
 569                 }
 570             }
 571         }
 572     }
 573
 574   /* Centre the moments */
 575   for ( j = 0 ; j < cov->dim - 1; ++j)
 576     {
 577       for (i = j + 1 ; i < cov->dim; ++i)
 578         {
 579           double *x = &cov->cm [cm_idx (cov, i, j)];
 580
 581           *x /= gsl_matrix_get (cov->moments[0], i, j);
 582
 583           *x -=
 584             gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)
 585             *
 586             gsl_matrix_get (cov->moments[MOMENT_MEAN], j, i);
 587         }
 588     }
 589
 590   return cm_to_gsl (cov);
 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   Covariance computed without dividing by the sample size.
 619  */
 620 static const gsl_matrix *
 621 covariance_calculate_double_pass_unnormalized (struct covariance *cov)
 622 {
 623   size_t i, j;
 624   for (i = 0 ; i < cov->dim; ++i)
 625     {
 626       for (j = 0 ; j < cov->dim; ++j)
 627         {
 628           int idx;
 629           double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
 630
 631           idx = cm_idx (cov, i, j);
 632           if ( idx >= 0)
 633             {
 634               x = &cov->cm [idx];
 635             }
 636         }
 637     }
 638
 639   return  cm_to_gsl (cov);
 640 }
 641
 642 static const gsl_matrix *
 643 covariance_calculate_single_pass_unnormalized (struct covariance *cov)
 644 {
 645   size_t i, j;
 646   size_t m;
 647
 648   for (i = 0 ; i < cov->dim; ++i)
 649     {
 650       for (j = 0 ; j < cov->dim; ++j)
 651         {
 652           double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
 653           *x -= pow2 (gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
 654             / gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 655         }
 656     }
 657   for ( j = 0 ; j < cov->dim - 1; ++j)
 658     {
 659       for (i = j + 1 ; i < cov->dim; ++i)
 660         {
 661           double *x = &cov->cm [cm_idx (cov, i, j)];
 662
 663           *x -=
 664             gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)
 665             *
 666             gsl_matrix_get (cov->moments[MOMENT_MEAN], j, i)
 667           / gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 668         }
 669     }
 670
 671   return cm_to_gsl (cov);
 672 }
 673
 674
 675 /*
 676    Return a pointer to gsl_matrix containing the pairwise covariances.
 677    The matrix remains owned by the COV object, and must not be freed.
 678    Call this function only after all data have been accumulated.
 679 */
 680 const gsl_matrix *
 681 covariance_calculate_unnormalized (struct covariance *cov)
 682 {
 683   assert ( cov->state > 0 );
 684
 685   switch (cov->passes)
 686     {
 687     case 1:
 688       return covariance_calculate_single_pass_unnormalized (cov);
 689       break;
 690     case 2:
 691       return covariance_calculate_double_pass_unnormalized (cov);
 692       break;
 693     default:
 694       NOT_REACHED ();
 695     }
 696 }
 697
 698
 699
 700 /* Destroy the COV object */
 701 void
 702 covariance_destroy (struct covariance *cov)
 703 {
 704   size_t i;
 705
 706   categoricals_destroy (cov->categoricals);
 707
 708   for (i = 0; i < n_MOMENTS; ++i)
 709     gsl_matrix_free (cov->moments[i]);
 710
 711   free (cov->moments);
 712   free (cov->cm);
 713   free (cov);
 714 }