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);
 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 static void
 236 dump_matrix (const gsl_matrix *m)
 237 {
 238   size_t i, j;
 239
 240   for (i = 0 ; i < m->size1; ++i)
 241     {
 242       for (j = 0 ; j < m->size2; ++j)
 243         printf ("%02f ", gsl_matrix_get (m, i, j));
 244       printf ("\n");
 245     }
 246 }
 247
 248 /* Call this function for every case in the data set */
 249 void
 250 covariance_accumulate_pass1 (struct covariance *cov, const struct ccase *c)
 251 {
 252   size_t i, j, m;
 253   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
 254
 255   assert (cov->passes == 2);
 256   if (!cov->pass_one_first_case_seen)
 257     {
 258       assert (cov->state == 0);
 259       cov->state = 1;
 260     }
 261
 262   categoricals_update (cov->categoricals, c);
 263
 264   for (i = 0 ; i < cov->dim; ++i)
 265     {
 266       const union value *val1 = case_data (c, cov->vars[i]);
 267
 268       if ( var_is_value_missing (cov->vars[i], val1, cov->exclude))
 269         continue;
 270
 271       for (j = 0 ; j < cov->dim; ++j)
 272         {
 273           double pwr = 1.0;
 274           const union value *val2 = case_data (c, cov->vars[j]);
 275
 276           if ( var_is_value_missing (cov->vars[j], val2, cov->exclude))
 277             continue;
 278
 279           for (m = 0 ; m <= MOMENT_MEAN; ++m)
 280             {
 281               double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 282
 283               *x += pwr * weight;
 284               pwr *= val1->f;
 285             }
 286         }
 287     }
 288
 289   cov->pass_one_first_case_seen = true;
 290 }
 291
 292
 293 /* Call this function for every case in the data set */
 294 void
 295 covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
 296 {
 297   size_t i, j;
 298   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
 299
 300   assert (cov->passes == 2);
 301   assert (cov->state >= 1);
 302
 303   if (! cov->pass_two_first_case_seen)
 304     {
 305       assert (cov->state == 1);
 306       cov->state = 2;
 307
 308       cov->dim = cov->n_vars + categoricals_total (cov->categoricals);
 309       cov->n_cm = (cov->dim * (cov->dim - 1)  ) / 2;
 310       cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm);
 311
 312       /* Grow the moment matrices so that they're large enough to accommodate the
 313          categorical elements */
 314       for (i = 0; i < n_MOMENTS; ++i)
 315         {
 316           cov->moments[i] = resize_matrix (cov->moments[i], cov->dim);
 317         }
 318
 319       /* Divide the means by the number of samples */
 320       for (i = 0; i < cov->n_vars; ++i)
 321         {
 322           for (j = 0; j < cov->n_vars; ++j)
 323             {
 324               double *x = gsl_matrix_ptr (cov->moments[MOMENT_MEAN], i, j);
 325               *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 326             }
 327         }
 328     }
 329
 330   for (i = 0 ; i < cov->dim; ++i)
 331     {
 332       const union value *val1 = case_data (c, cov->vars[i]);
 333
 334       if ( var_is_value_missing (cov->vars[i], val1, cov->exclude))
 335         continue;
 336
 337       for (j = 0 ; j < cov->dim; ++j)
 338         {
 339           int idx;
 340           double ss ;
 341           const union value *val2 = case_data (c, cov->vars[j]);
 342
 343           const double s = pow2 (val1->f - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)) * weight;
 344
 345           if ( var_is_value_missing (cov->vars[j], val2, cov->exclude))
 346             continue;
 347
 348           {
 349             double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
 350             *x += s;
 351           }
 352
 353           ss =
 354             (val1->f - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
 355             *
 356             (val2->f - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
 357             * weight
 358             ;
 359
 360           idx = cm_idx (cov, i, j);
 361           if (idx >= 0)
 362             {
 363               cov->cm [idx] += ss;
 364             }
 365
 366         }
 367     }
 368
 369   cov->pass_two_first_case_seen = true;
 370 }
 371
 372
 373 /* Call this function for every case in the data set.
 374    After all cases have been passed, call covariance_calculate
 375  */
 376 void
 377 covariance_accumulate (struct covariance *cov, const struct ccase *c)
 378 {
 379   size_t i, j, m;
 380   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
 381
 382   assert (cov->passes == 1);
 383
 384   if ( !cov->pass_one_first_case_seen)
 385     {
 386       assert ( cov->state == 0);
 387       cov->state = 1;
 388     }
 389
 390   for (i = 0 ; i < cov->dim; ++i)
 391     {
 392       const union value *val1 = case_data (c, cov->vars[i]);
 393
 394       if ( var_is_value_missing (cov->vars[i], val1, cov->exclude))
 395         continue;
 396
 397       for (j = 0 ; j < cov->dim; ++j)
 398         {
 399           double pwr = 1.0;
 400           int idx;
 401           const union value *val2 = case_data (c, cov->vars[j]);
 402
 403           if ( var_is_value_missing (cov->vars[j], val2, cov->exclude))
 404             continue;
 405
 406           idx = cm_idx (cov, i, j);
 407           if (idx >= 0)
 408             {
 409               cov->cm [idx] += val1->f * val2->f * weight;
 410             }
 411
 412           for (m = 0 ; m < n_MOMENTS; ++m)
 413             {
 414               double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 415
 416               *x += pwr * weight;
 417               pwr *= val1->f;
 418             }
 419         }
 420     }
 421
 422   cov->pass_one_first_case_seen = true;
 423 }
 424
 425
 426 /*
 427    Allocate and return a gsl_matrix containing the covariances of the
 428    data.
 429 */
 430 static gsl_matrix *
 431 cm_to_gsl (struct covariance *cov)
 432 {
 433   int i, j;
 434   gsl_matrix *m = gsl_matrix_calloc (cov->dim, cov->dim);
 435
 436   /* Copy the non-diagonal elements from cov->cm */
 437   for ( j = 0 ; j < cov->dim - 1; ++j)
 438     {
 439       for (i = j+1 ; i < cov->dim; ++i)
 440         {
 441           double x = cov->cm [cm_idx (cov, i, j)];
 442           gsl_matrix_set (m, i, j, x);
 443           gsl_matrix_set (m, j, i, x);
 444         }
 445     }
 446
 447   /* Copy the diagonal elements from cov->moments[2] */
 448   for (j = 0 ; j < cov->dim ; ++j)
 449     {
 450       double sigma = gsl_matrix_get (cov->moments[2], j, j);
 451       gsl_matrix_set (m, j, j, sigma);
 452     }
 453
 454   return m;
 455 }
 456
 457
 458 static const gsl_matrix *
 459 covariance_calculate_double_pass (struct covariance *cov)
 460 {
 461   size_t i, j;
 462   for (i = 0 ; i < cov->dim; ++i)
 463     {
 464       for (j = 0 ; j < cov->dim; ++j)
 465         {
 466           int idx;
 467           double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
 468           *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 469
 470           idx = cm_idx (cov, i, j);
 471           if ( idx >= 0)
 472             {
 473               x = &cov->cm [idx];
 474               *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 475             }
 476         }
 477     }
 478
 479   return  cm_to_gsl (cov);
 480 }
 481
 482 static const gsl_matrix *
 483 covariance_calculate_single_pass (struct covariance *cov)
 484 {
 485   size_t i, j;
 486   size_t m;
 487
 488   for (m = 0; m < n_MOMENTS; ++m)
 489     {
 490       /* Divide the moments by the number of samples */
 491       if ( m > 0)
 492         {
 493           for (i = 0 ; i < cov->dim; ++i)
 494             {
 495               for (j = 0 ; j < cov->dim; ++j)
 496                 {
 497                   double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 498                   *x /= gsl_matrix_get (cov->moments[0], i, j);
 499
 500                   if ( m == MOMENT_VARIANCE)
 501                     *x -= pow2 (gsl_matrix_get (cov->moments[1], i, j));
 502                 }
 503             }
 504         }
 505     }
 506
 507   /* Centre the moments */
 508   for ( j = 0 ; j < cov->dim - 1; ++j)
 509     {
 510       for (i = j + 1 ; i < cov->dim; ++i)
 511         {
 512           double *x = &cov->cm [cm_idx (cov, i, j)];
 513
 514           *x /= gsl_matrix_get (cov->moments[0], i, j);
 515
 516           *x -=
 517             gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)
 518             *
 519             gsl_matrix_get (cov->moments[MOMENT_MEAN], j, i);
 520         }
 521     }
 522
 523   return cm_to_gsl (cov);
 524 }
 525
 526
 527
 528 /*
 529    Return a pointer to gsl_matrix containing the pairwise covariances.
 530    The matrix remains owned by the COV object, and must not be freed.
 531    Call this function only after all data have been accumulated.
 532 */
 533 const gsl_matrix *
 534 covariance_calculate (struct covariance *cov)
 535 {
 536   assert ( cov->state > 0 );
 537
 538   switch (cov->passes)
 539     {
 540     case 1:
 541       return covariance_calculate_single_pass (cov);
 542       break;
 543     case 2:
 544       return covariance_calculate_double_pass (cov);
 545       break;
 546     default:
 547       NOT_REACHED ();
 548     }
 549 }
 550
 551
 552
 553
 554 /* Destroy the COV object */
 555 void
 556 covariance_destroy (struct covariance *cov)
 557 {
 558   size_t i;
 559   free (cov->vars);
 560   categoricals_destroy (cov->categoricals);
 561
 562   for (i = 0; i < n_MOMENTS; ++i)
 563     gsl_matrix_free (cov->moments[i]);
 564
 565   free (cov->moments);
 566   free (cov->cm);
 567   free (cov);
 568 }