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