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