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
  28 #define n_MOMENTS (MOMENT_VARIANCE + 1)
  29
  30
  31 struct covariance
  32 {
  33   /* The variables for which the covariance matrix is to be calculated. */
  34   size_t n_vars;
  35   const struct variable **vars;
  36
  37   /* Categorical variables. */
  38   size_t n_catvars;
  39   const struct variable **catvars;
  40
  41   /* Array containing number of categories per categorical variable. */
  42   size_t *n_categories;
  43
  44   /* Dimension of the covariance matrix. */
  45   size_t dim;
  46
  47   /* The weight variable (or NULL if none) */
  48   const struct variable *wv;
  49
  50   /* A set of matrices containing the 0th, 1st and 2nd moments */
  51   gsl_matrix **moments;
  52
  53   /* The class of missing values to exclude */
  54   enum mv_class exclude;
  55
  56   /* An array of doubles representing the covariance matrix.
  57      Only the top triangle is included, and no diagonals */
  58   double *cm;
  59   int n_cm;
  60
  61   /* 1 for single pass algorithm;
  62      2 for double pass algorithm
  63   */
  64   short passes;
  65
  66   /*
  67     0 : No pass has  been made
  68     1 : First pass has been started
  69     2 : Second pass has been
  70
  71     IE: How many passes have been (partially) made. */
  72   short state;
  73
  74   /* Flags indicating that the first case has been seen */
  75   bool pass_one_first_case_seen;
  76   bool pass_two_first_case_seen;
  77 };
  78
  79
  80
  81 /* Return a matrix containing the M th moments.
  82    The matrix is of size  NxN where N is the number of variables.
  83    Each row represents the moments of a variable.
  84    In the absence of missing values, the columns of this matrix will
  85    be identical.  If missing values are involved, then element (i,j)
  86    is the moment of the i th variable, when paired with the j th variable.
  87  */
  88 const gsl_matrix *
  89 covariance_moments (const struct covariance *cov, int m)
  90 {
  91   return cov->moments[m];
  92 }
  93
  94
  95
  96 /* Create a covariance struct.
  97  */
  98 struct covariance *
  99 covariance_create (size_t n_vars, const struct variable **vars,
 100                    const struct variable *weight, enum mv_class exclude,
 101                    short passes)
 102 {
 103   size_t i;
 104   struct covariance *cov = xmalloc (sizeof *cov);
 105   assert (passes == 1 || passes == 2);
 106   cov->passes = passes;
 107   cov->state = 0;
 108   cov->pass_one_first_case_seen = cov->pass_two_first_case_seen = false;
 109
 110   cov->vars = xmalloc (sizeof *cov->vars * n_vars);
 111
 112   cov->wv = weight;
 113   cov->n_vars = n_vars;
 114   cov->dim = n_vars;
 115
 116   for (i = 0; i < n_vars; ++i)
 117     cov->vars[i] = vars[i];
 118
 119   cov->moments = xmalloc (sizeof *cov->moments * n_MOMENTS);
 120
 121   for (i = 0; i < n_MOMENTS; ++i)
 122     cov->moments[i] = gsl_matrix_calloc (n_vars, n_vars);
 123
 124   cov->exclude = exclude;
 125
 126   cov->n_cm = (n_vars * (n_vars - 1)  ) / 2;
 127
 128   cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm);
 129
 130   return cov;
 131 }
 132
 133 /*
 134   Create a covariance struct for a two-pass algorithm. If categorical
 135   variables are involed, the dimension cannot be know until after the
 136   first data pass, so the actual covariances will not be allocated
 137   until then.
 138  */
 139 struct covariance *
 140 covariance_2pass_create (size_t n_vars, const struct variable **vars,
 141                          size_t n_catvars, const struct variable **catvars,
 142                          const struct variable *weight, enum mv_class exclude)
 143 {
 144   size_t i;
 145   struct covariance *cov = xmalloc (sizeof *cov);
 146   cov->vars = xmalloc (sizeof *cov->vars * n_vars);
 147   cov->catvars = xnmalloc (n_catvars, sizeof (*cov->catvars));
 148   cov->n_categories = xnmalloc (n_catvars, sizeof (cov->n_categories));
 149
 150   cov->wv = weight;
 151   cov->n_vars = n_vars;
 152   cov->n_catvars = n_catvars;
 153
 154   for (i = 0; i < n_vars; ++i)
 155     cov->vars[i] = vars[i];
 156
 157   for (i = 0; i < n_catvars; i++)
 158     {
 159       cov->catvars[i] = catvars[i];
 160       cov->n_categories[i] = 0;
 161     }
 162
 163   cov->moments = xmalloc (sizeof *cov->moments * n_MOMENTS);
 164
 165   cov->exclude = exclude;
 166
 167   return cov;
 168 }
 169
 170 /* Return an integer, which can be used to index
 171    into COV->cm, to obtain the I, J th element
 172    of the covariance matrix.  If COV->cm does not
 173    contain that element, then a negative value
 174    will be returned.
 175 */
 176 static int
 177 cm_idx (const struct covariance *cov, int i, int j)
 178 {
 179   int as;
 180   const int n2j = cov->n_vars - 2 - j;
 181   const int nj = cov->n_vars - 2 ;
 182
 183   assert (i >= 0);
 184   assert (j < cov->n_vars);
 185
 186   if ( i == 0)
 187     return -1;
 188
 189   if (j >= cov->n_vars - 1)
 190     return -1;
 191
 192   if ( i <= j)
 193     return -1 ;
 194
 195   as = nj * (nj + 1) ;
 196   as -= n2j * (n2j + 1) ;
 197   as /= 2;
 198
 199   return i - 1 + as;
 200 }
 201
 202 static void
 203 dump_matrix (const gsl_matrix *m)
 204 {
 205   size_t i, j;
 206
 207   for (i = 0 ; i < m->size1; ++i)
 208     {
 209       for (j = 0 ; j < m->size2; ++j)
 210         printf ("%02f ", gsl_matrix_get (m, i, j));
 211       printf ("\n");
 212     }
 213 }
 214
 215 /* Call this function for every case in the data set */
 216 void
 217 covariance_accumulate_pass1 (struct covariance *cov, const struct ccase *c)
 218 {
 219   size_t i, j, m;
 220   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
 221
 222   assert (cov->passes == 2);
 223   if (!cov->pass_one_first_case_seen)
 224     {
 225       assert (cov->state == 0);
 226       cov->state = 1;
 227     }
 228
 229   for (i = 0 ; i < cov->n_vars; ++i)
 230     {
 231       const union value *val1 = case_data (c, cov->vars[i]);
 232
 233       if ( var_is_value_missing (cov->vars[i], val1, cov->exclude))
 234         continue;
 235
 236       for (j = 0 ; j < cov->n_vars; ++j)
 237         {
 238           double pwr = 1.0;
 239           const union value *val2 = case_data (c, cov->vars[j]);
 240
 241           if ( var_is_value_missing (cov->vars[j], val2, cov->exclude))
 242             continue;
 243
 244           for (m = 0 ; m <= MOMENT_MEAN; ++m)
 245             {
 246               double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 247
 248               *x += pwr * weight;
 249               pwr *= val1->f;
 250             }
 251         }
 252     }
 253
 254   cov->pass_one_first_case_seen = true;
 255 }
 256
 257
 258 /* Call this function for every case in the data set */
 259 void
 260 covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
 261 {
 262   size_t i, j;
 263   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
 264
 265   assert (cov->passes == 2);
 266   assert (cov->state >= 1);
 267
 268   if (! cov->pass_two_first_case_seen)
 269     {
 270       assert (cov->state == 1);
 271       cov->state = 2;
 272
 273       /* Divide the means by the number of samples */
 274       for (i = 0; i < cov->n_vars; ++i)
 275         {
 276           for (j = 0; j < cov->n_vars; ++j)
 277             {
 278               double *x = gsl_matrix_ptr (cov->moments[MOMENT_MEAN], i, j);
 279               *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 280             }
 281         }
 282     }
 283
 284   for (i = 0 ; i < cov->n_vars; ++i)
 285     {
 286       const union value *val1 = case_data (c, cov->vars[i]);
 287
 288       if ( var_is_value_missing (cov->vars[i], val1, cov->exclude))
 289         continue;
 290
 291       for (j = 0 ; j < cov->n_vars; ++j)
 292         {
 293           int idx;
 294           double ss ;
 295           const union value *val2 = case_data (c, cov->vars[j]);
 296
 297           const double s = pow2 (val1->f - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)) * weight;
 298
 299           if ( var_is_value_missing (cov->vars[j], val2, cov->exclude))
 300             continue;
 301
 302           {
 303             double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
 304             *x += s;
 305           }
 306
 307           ss =
 308             (val1->f - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
 309             *
 310             (val2->f - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
 311             * weight
 312             ;
 313
 314           idx = cm_idx (cov, i, j);
 315           if (idx >= 0)
 316             {
 317               cov->cm [idx] += ss;
 318             }
 319
 320         }
 321     }
 322
 323   cov->pass_two_first_case_seen = true;
 324 }
 325
 326
 327 /* Call this function for every case in the data set.
 328    After all cases have been passed, call covariance_calculate
 329  */
 330 void
 331 covariance_accumulate (struct covariance *cov, const struct ccase *c)
 332 {
 333   size_t i, j, m;
 334   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
 335
 336   assert (cov->passes == 1);
 337
 338   if ( !cov->pass_one_first_case_seen)
 339     {
 340       assert ( cov->state == 0);
 341       cov->state = 1;
 342     }
 343
 344   for (i = 0 ; i < cov->n_vars; ++i)
 345     {
 346       const union value *val1 = case_data (c, cov->vars[i]);
 347
 348       if ( var_is_value_missing (cov->vars[i], val1, cov->exclude))
 349         continue;
 350
 351       for (j = 0 ; j < cov->n_vars; ++j)
 352         {
 353           double pwr = 1.0;
 354           int idx;
 355           const union value *val2 = case_data (c, cov->vars[j]);
 356
 357           if ( var_is_value_missing (cov->vars[j], val2, cov->exclude))
 358             continue;
 359
 360           idx = cm_idx (cov, i, j);
 361           if (idx >= 0)
 362             {
 363               cov->cm [idx] += val1->f * val2->f * weight;
 364             }
 365
 366           for (m = 0 ; m < n_MOMENTS; ++m)
 367             {
 368               double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 369
 370               *x += pwr * weight;
 371               pwr *= val1->f;
 372             }
 373         }
 374     }
 375
 376   cov->pass_one_first_case_seen = true;
 377 }
 378
 379
 380 /*
 381    Allocate and return a gsl_matrix containing the covariances of the
 382    data.
 383 */
 384 static gsl_matrix *
 385 cm_to_gsl (struct covariance *cov)
 386 {
 387   int i, j;
 388   gsl_matrix *m = gsl_matrix_calloc (cov->n_vars, cov->n_vars);
 389
 390   /* Copy the non-diagonal elements from cov->cm */
 391   for ( j = 0 ; j < cov->n_vars - 1; ++j)
 392     {
 393       for (i = j+1 ; i < cov->n_vars; ++i)
 394         {
 395           double x = cov->cm [cm_idx (cov, i, j)];
 396           gsl_matrix_set (m, i, j, x);
 397           gsl_matrix_set (m, j, i, x);
 398         }
 399     }
 400
 401   /* Copy the diagonal elements from cov->moments[2] */
 402   for (j = 0 ; j < cov->n_vars ; ++j)
 403     {
 404       double sigma = gsl_matrix_get (cov->moments[2], j, j);
 405       gsl_matrix_set (m, j, j, sigma);
 406     }
 407
 408   return m;
 409 }
 410
 411
 412 static const gsl_matrix *
 413 covariance_calculate_double_pass (struct covariance *cov)
 414 {
 415   size_t i, j;
 416   for (i = 0 ; i < cov->n_vars; ++i)
 417     {
 418       for (j = 0 ; j < cov->n_vars; ++j)
 419         {
 420           int idx;
 421           double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
 422           *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 423
 424           idx = cm_idx (cov, i, j);
 425           if ( idx >= 0)
 426             {
 427               x = &cov->cm [idx];
 428               *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
 429             }
 430         }
 431     }
 432
 433   return  cm_to_gsl (cov);
 434 }
 435
 436 static const gsl_matrix *
 437 covariance_calculate_single_pass (struct covariance *cov)
 438 {
 439   size_t i, j;
 440   size_t m;
 441
 442   for (m = 0; m < n_MOMENTS; ++m)
 443     {
 444       /* Divide the moments by the number of samples */
 445       if ( m > 0)
 446         {
 447           for (i = 0 ; i < cov->n_vars; ++i)
 448             {
 449               for (j = 0 ; j < cov->n_vars; ++j)
 450                 {
 451                   double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 452                   *x /= gsl_matrix_get (cov->moments[0], i, j);
 453
 454                   if ( m == MOMENT_VARIANCE)
 455                     *x -= pow2 (gsl_matrix_get (cov->moments[1], i, j));
 456                 }
 457             }
 458         }
 459     }
 460
 461   /* Centre the moments */
 462   for ( j = 0 ; j < cov->n_vars - 1; ++j)
 463     {
 464       for (i = j + 1 ; i < cov->n_vars; ++i)
 465         {
 466           double *x = &cov->cm [cm_idx (cov, i, j)];
 467
 468           *x /= gsl_matrix_get (cov->moments[0], i, j);
 469
 470           *x -=
 471             gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)
 472             *
 473             gsl_matrix_get (cov->moments[MOMENT_MEAN], j, i);
 474         }
 475     }
 476
 477   return cm_to_gsl (cov);
 478 }
 479
 480
 481
 482 /*
 483    Return a pointer to gsl_matrix containing the pairwise covariances.
 484    The matrix remains owned by the COV object, and must not be freed.
 485    Call this function only after all data have been accumulated.
 486 */
 487 const gsl_matrix *
 488 covariance_calculate (struct covariance *cov)
 489 {
 490   assert ( cov->state > 0 );
 491
 492   switch (cov->passes)
 493     {
 494     case 1:
 495       return covariance_calculate_single_pass (cov);
 496       break;
 497     case 2:
 498       return covariance_calculate_double_pass (cov);
 499       break;
 500     default:
 501       NOT_REACHED ();
 502     }
 503 }
 504
 505
 506
 507
 508 /* Destroy the COV object */
 509 void
 510 covariance_destroy (struct covariance *cov)
 511 {
 512   size_t i;
 513   free (cov->vars);
 514
 515   for (i = 0; i < n_MOMENTS; ++i)
 516     gsl_matrix_free (cov->moments[i]);
 517
 518   free (cov->moments);
 519   free (cov->cm);
 520   free (cov);
 521 }