src/math/covariance.c

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