pintos-os.org Git - pspp/blob - 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 *const *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 *const *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 *const *vars,
 174                          struct categoricals *cats,
 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 = cats;
 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 #if 0
 271 void
 272 dump_matrix (const gsl_matrix *m)
 273 {
 274   size_t i, j;
 275
 276   for (i = 0 ; i < m->size1; ++i)
 277     {
 278       for (j = 0 ; j < m->size2; ++j)
 279         printf ("%02f ", gsl_matrix_get (m, i, j));
 280       printf ("\n");
 281     }
 282 }
 283 #endif
 284
 285 /* Call this function for every case in the data set */
 286 void
 287 covariance_accumulate_pass1 (struct covariance *cov, const struct ccase *c)
 288 {
 289   size_t i, j, m;
 290   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
 291
 292   assert (cov->passes == 2);
 293   if (!cov->pass_one_first_case_seen)
 294     {
 295       assert (cov->state == 0);
 296       cov->state = 1;
 297     }
 298
 299   if (cov->categoricals)
 300     categoricals_update (cov->categoricals, c);
 301
 302   for (i = 0 ; i < cov->dim; ++i)
 303     {
 304       double v1 = get_val (cov, i, c);
 305
 306       if ( is_missing (cov, i, c))
 307         continue;
 308
 309       for (j = 0 ; j < cov->dim; ++j)
 310         {
 311           double pwr = 1.0;
 312
 313           if ( is_missing (cov, j, c))
 314             continue;
 315
 316           for (m = 0 ; m <= MOMENT_MEAN; ++m)
 317             {
 318               double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 319
 320               *x += pwr * weight;
 321               pwr *= v1;
 322             }
 323         }
 324     }
 325
 326   cov->pass_one_first_case_seen = true;
 327 }
 328
 329
 330 /* Call this function for every case in the data set */
 331 void
 332 covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
 333 {
 334   size_t i, j;
 335   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
 336
 337   assert (cov->passes == 2);
 338   assert (cov->state >= 1);
 339
 340   if (! cov->pass_two_first_case_seen)
 341     {
 342       size_t m;
 343       assert (cov->state == 1);
 344       cov->state = 2;
 345
 346       cov->dim = cov->n_vars;
 347
 348       if (cov->categoricals)
 349         cov->dim += categoricals_total (cov->categoricals)
 350           - categoricals_get_n_variables (cov->categoricals);
 351
 352       cov->n_cm = (cov->dim * (cov->dim - 1)  ) / 2;
 353       cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm);
 354
 355       /* Grow the moment matrices so that they're large enough to accommodate the
 356          categorical elements */
 357       for (i = 0; i < n_MOMENTS; ++i)
 358         {
 359           cov->moments[i] = resize_matrix (cov->moments[i], cov->dim);
 360         }
 361
 362       if (cov->categoricals)
 363         categoricals_done (cov->categoricals);
 364
 365       /* Populate the moments matrices with the categorical value elements */
 366       for (i = cov->n_vars; i < cov->dim; ++i)
 367         {
 368           for (j = 0 ; j < cov->dim ; ++j) /* FIXME: This is WRONG !!! */
 369             {
 370               double w = categoricals_get_weight_by_subscript (cov->categoricals, i - cov->n_vars);
 371
 372               gsl_matrix_set (cov->moments[MOMENT_NONE], i, j, w);
 373
 374               w = categoricals_get_sum_by_subscript (cov->categoricals, i - cov->n_vars);
 375
 376               gsl_matrix_set (cov->moments[MOMENT_MEAN], i, j, w);
 377             }
 378         }
 379
 380       /* FIXME: This is WRONG!!  It must be fixed to properly handle missing values.  For
 381        now it assumes there are none */
 382       for (m = 0 ; m < n_MOMENTS; ++m)
 383         {
 384           for (i = 0 ; i < cov->dim ; ++i)
 385             {
 386               double x = gsl_matrix_get (cov->moments[m], i, cov->n_vars -1);
 387               for (j = cov->n_vars; j < cov->dim; ++j)
 388                 {
 389                   gsl_matrix_set (cov->moments[m], i, j, x);
 390                 }
 391             }
 392         }
 393
 394       /* Divide the means by the number of samples */
 395       for (i = 0; i < cov->dim; ++i)
 396         {
 397           for (j = 0; j < cov->dim; ++j)
 398             {
 399               double *x = gsl_matrix_ptr (cov->moments[MOMENT_MEAN], i, j);
 400               *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 401             }
 402         }
 403     }
 404
 405   for (i = 0 ; i < cov->dim; ++i)
 406     {
 407       double v1 = get_val (cov, i, c);
 408
 409       if ( is_missing (cov, i, c))
 410         continue;
 411
 412       for (j = 0 ; j < cov->dim; ++j)
 413         {
 414           int idx;
 415           double ss ;
 416           double v2 = get_val (cov, j, c);
 417
 418           const double s = pow2 (v1 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)) * weight;
 419
 420           if ( is_missing (cov, j, c))
 421             continue;
 422
 423           {
 424             double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
 425             *x += s;
 426           }
 427
 428           ss =
 429             (v1 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
 430             *
 431             (v2 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
 432             * weight
 433             ;
 434
 435           idx = cm_idx (cov, i, j);
 436           if (idx >= 0)
 437             {
 438               cov->cm [idx] += ss;
 439             }
 440         }
 441     }
 442
 443   cov->pass_two_first_case_seen = true;
 444 }
 445
 446
 447 /* Call this function for every case in the data set.
 448    After all cases have been passed, call covariance_calculate
 449  */
 450 void
 451 covariance_accumulate (struct covariance *cov, const struct ccase *c)
 452 {
 453   size_t i, j, m;
 454   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
 455
 456   assert (cov->passes == 1);
 457
 458   if ( !cov->pass_one_first_case_seen)
 459     {
 460       assert ( cov->state == 0);
 461       cov->state = 1;
 462     }
 463
 464   for (i = 0 ; i < cov->dim; ++i)
 465     {
 466       const union value *val1 = case_data (c, cov->vars[i]);
 467
 468       if ( is_missing (cov, i, c))
 469         continue;
 470
 471       for (j = 0 ; j < cov->dim; ++j)
 472         {
 473           double pwr = 1.0;
 474           int idx;
 475           const union value *val2 = case_data (c, cov->vars[j]);
 476
 477           if ( is_missing (cov, j, c))
 478             continue;
 479
 480           idx = cm_idx (cov, i, j);
 481           if (idx >= 0)
 482             {
 483               cov->cm [idx] += val1->f * val2->f * weight;
 484             }
 485
 486           for (m = 0 ; m < n_MOMENTS; ++m)
 487             {
 488               double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 489
 490               *x += pwr * weight;
 491               pwr *= val1->f;
 492             }
 493         }
 494     }
 495
 496   cov->pass_one_first_case_seen = true;
 497 }
 498
 499
 500 /*
 501    Allocate and return a gsl_matrix containing the covariances of the
 502    data.
 503 */
 504 static gsl_matrix *
 505 cm_to_gsl (struct covariance *cov)
 506 {
 507   int i, j;
 508   gsl_matrix *m = gsl_matrix_calloc (cov->dim, cov->dim);
 509
 510   /* Copy the non-diagonal elements from cov->cm */
 511   for ( j = 0 ; j < cov->dim - 1; ++j)
 512     {
 513       for (i = j+1 ; i < cov->dim; ++i)
 514         {
 515           double x = cov->cm [cm_idx (cov, i, j)];
 516           gsl_matrix_set (m, i, j, x);
 517           gsl_matrix_set (m, j, i, x);
 518         }
 519     }
 520
 521   /* Copy the diagonal elements from cov->moments[2] */
 522   for (j = 0 ; j < cov->dim ; ++j)
 523     {
 524       double sigma = gsl_matrix_get (cov->moments[2], j, j);
 525       gsl_matrix_set (m, j, j, sigma);
 526     }
 527
 528   return m;
 529 }
 530
 531
 532 static const gsl_matrix *
 533 covariance_calculate_double_pass (struct covariance *cov)
 534 {
 535   size_t i, j;
 536   for (i = 0 ; i < cov->dim; ++i)
 537     {
 538       for (j = 0 ; j < cov->dim; ++j)
 539         {
 540           int idx;
 541           double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
 542           *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 543
 544           idx = cm_idx (cov, i, j);
 545           if ( idx >= 0)
 546             {
 547               x = &cov->cm [idx];
 548               *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 549             }
 550         }
 551     }
 552
 553   return  cm_to_gsl (cov);
 554 }
 555
 556 static const gsl_matrix *
 557 covariance_calculate_single_pass (struct covariance *cov)
 558 {
 559   size_t i, j;
 560   size_t m;
 561
 562   for (m = 0; m < n_MOMENTS; ++m)
 563     {
 564       /* Divide the moments by the number of samples */
 565       if ( m > 0)
 566         {
 567           for (i = 0 ; i < cov->dim; ++i)
 568             {
 569               for (j = 0 ; j < cov->dim; ++j)
 570                 {
 571                   double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 572                   *x /= gsl_matrix_get (cov->moments[0], i, j);
 573
 574                   if ( m == MOMENT_VARIANCE)
 575                     *x -= pow2 (gsl_matrix_get (cov->moments[1], i, j));
 576                 }
 577             }
 578         }
 579     }
 580
 581   /* Centre the moments */
 582   for ( j = 0 ; j < cov->dim - 1; ++j)
 583     {
 584       for (i = j + 1 ; i < cov->dim; ++i)
 585         {
 586           double *x = &cov->cm [cm_idx (cov, i, j)];
 587
 588           *x /= gsl_matrix_get (cov->moments[0], i, j);
 589
 590           *x -=
 591             gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)
 592             *
 593             gsl_matrix_get (cov->moments[MOMENT_MEAN], j, i);
 594         }
 595     }
 596
 597   return cm_to_gsl (cov);
 598 }
 599
 600
 601 /*
 602    Return a pointer to gsl_matrix containing the pairwise covariances.
 603    The matrix remains owned by the COV object, and must not be freed.
 604    Call this function only after all data have been accumulated.
 605 */
 606 const gsl_matrix *
 607 covariance_calculate (struct covariance *cov)
 608 {
 609   if ( cov->state <= 0 )
 610     return NULL;
 611
 612   switch (cov->passes)
 613     {
 614     case 1:
 615       return covariance_calculate_single_pass (cov);
 616       break;
 617     case 2:
 618       return covariance_calculate_double_pass (cov);
 619       break;
 620     default:
 621       NOT_REACHED ();
 622     }
 623 }
 624
 625 /*
 626   Covariance computed without dividing by the sample size.
 627  */
 628 static const gsl_matrix *
 629 covariance_calculate_double_pass_unnormalized (struct covariance *cov)
 630 {
 631   size_t i, j;
 632   for (i = 0 ; i < cov->dim; ++i)
 633     {
 634       for (j = 0 ; j < cov->dim; ++j)
 635         {
 636           int idx;
 637           double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
 638
 639           idx = cm_idx (cov, i, j);
 640           if ( idx >= 0)
 641             {
 642               x = &cov->cm [idx];
 643             }
 644         }
 645     }
 646
 647   return  cm_to_gsl (cov);
 648 }
 649
 650 static const gsl_matrix *
 651 covariance_calculate_single_pass_unnormalized (struct covariance *cov)
 652 {
 653   size_t i, j;
 654
 655   for (i = 0 ; i < cov->dim; ++i)
 656     {
 657       for (j = 0 ; j < cov->dim; ++j)
 658         {
 659           double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
 660           *x -= pow2 (gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
 661             / gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 662         }
 663     }
 664   for ( j = 0 ; j < cov->dim - 1; ++j)
 665     {
 666       for (i = j + 1 ; i < cov->dim; ++i)
 667         {
 668           double *x = &cov->cm [cm_idx (cov, i, j)];
 669
 670           *x -=
 671             gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)
 672             *
 673             gsl_matrix_get (cov->moments[MOMENT_MEAN], j, i)
 674           / gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 675         }
 676     }
 677
 678   return cm_to_gsl (cov);
 679 }
 680
 681
 682 /*
 683    Return a pointer to gsl_matrix containing the pairwise covariances.
 684    The matrix remains owned by the COV object, and must not be freed.
 685    Call this function only after all data have been accumulated.
 686 */
 687 const gsl_matrix *
 688 covariance_calculate_unnormalized (struct covariance *cov)
 689 {
 690   if ( cov->state <= 0 )
 691     return NULL;
 692
 693   switch (cov->passes)
 694     {
 695     case 1:
 696       return covariance_calculate_single_pass_unnormalized (cov);
 697       break;
 698     case 2:
 699       return covariance_calculate_double_pass_unnormalized (cov);
 700       break;
 701     default:
 702       NOT_REACHED ();
 703     }
 704 }
 705
 706 /* Function to access the categoricals used by COV
 707    The return value is owned by the COV
 708 */
 709 const struct categoricals *
 710 covariance_get_categoricals (const struct covariance *cov)
 711 {
 712   return cov->categoricals;
 713 }
 714
 715
 716 /* Destroy the COV object */
 717 void
 718 covariance_destroy (struct covariance *cov)
 719 {
 720   size_t i;
 721
 722   categoricals_destroy (cov->categoricals);
 723
 724   for (i = 0; i < n_MOMENTS; ++i)
 725     gsl_matrix_free (cov->moments[i]);
 726
 727   free (cov->moments);
 728   free (cov->cm);
 729   free (cov);
 730 }