pintos-os.org Git - pspp/blob - src/math/linreg.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2005, 2010, 2011, 2017 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/linreg.h"
  20
  21 #include <gsl/gsl_blas.h>
  22 #include <gsl/gsl_cblas.h>
  23 #include <gsl/gsl_errno.h>
  24 #include <gsl/gsl_fit.h>
  25 #include <gsl/gsl_linalg.h>
  26 #include <gsl/gsl_multifit.h>
  27
  28 #include "data/value.h"
  29 #include "data/variable.h"
  30 #include "linreg/sweep.h"
  31
  32 #include "gl/xalloc.h"
  33
  34 /*
  35   Find the least-squares estimate of b for the linear model:
  36
  37   Y = Xb + Z
  38
  39   where Y is an n-by-1 column vector, X is an n-by-p matrix of
  40   independent variables, b is a p-by-1 vector of regression coefficients,
  41   and Z is an n-by-1 normally-distributed random vector with independent
  42   identically distributed components with mean 0.
  43
  44   This estimate is found via the sweep operator or singular-value
  45   decomposition with gsl.
  46
  47
  48   References:
  49
  50   1. Matrix Computations, third edition. GH Golub and CF Van Loan.
  51   The Johns Hopkins University Press. 1996. ISBN 0-8018-5414-8.
  52
  53   2. Numerical Analysis for Statisticians. K Lange. Springer. 1999.
  54   ISBN 0-387-94979-8.
  55
  56   3. Numerical Linear Algebra for Applications in Statistics. JE Gentle.
  57   Springer. 1998. ISBN 0-387-98542-5.
  58 */
  59
  60 struct linreg
  61 {
  62   double n_obs;                 /* Number of observations. */
  63   int n_indeps;                 /* Number of independent variables. */
  64   int n_coeffs;                 /* The intercept is not considered a
  65                                    coefficient here. */
  66
  67   /*
  68     Pointers to the variables.
  69    */
  70   const struct variable *depvar;
  71   const struct variable **indep_vars;
  72
  73   double *coeff;
  74   double intercept;
  75   /*
  76      Means and standard deviations of the variables.
  77      If these pointers are null when pspp_linreg() is
  78      called, pspp_linreg() will compute their values.
  79
  80      Entry i of indep_means is the mean of independent
  81      variable i, whose observations are stored in the ith
  82      column of the design matrix.
  83    */
  84   double depvar_mean;
  85   gsl_vector *indep_means;
  86   gsl_vector *indep_std;
  87
  88   /*
  89      Sums of squares.
  90    */
  91   double ssm;                   /* Sums of squares for the overall model. */
  92   double sst;                   /* Sum of squares total. */
  93   double sse;                   /* Sum of squares error. */
  94   double mse;                   /* Mean squared error. This is just sse /
  95                                    dfe, but since it is the best unbiased
  96                                    estimate of the population variance, it
  97                                    has its own entry here. */
  98   /*
  99      Covariance matrix of the parameter estimates.
 100    */
 101   gsl_matrix *cov;
 102   /*
 103      Degrees of freedom.
 104    */
 105   double dft;
 106   double dfe;
 107   double dfm;
 108
 109   int dependent_column; /* Column containing the dependent variable. Defaults to last column. */
 110   int refcnt;
 111
 112   bool origin;
 113 };
 114
 115 const struct variable **
 116 linreg_get_vars (const struct linreg *c)
 117 {
 118   return c->indep_vars;
 119 }
 120
 121 /*
 122   Allocate a linreg and return a pointer to it. n is the number of
 123   cases, p is the number of independent variables.
 124  */
 125 struct linreg *
 126 linreg_alloc (const struct variable *depvar, const struct variable **indep_vars,
 127               double n, size_t p, bool origin)
 128 {
 129   struct linreg *c;
 130   size_t i;
 131
 132   c = xmalloc (sizeof (*c));
 133   c->depvar = depvar;
 134   c->indep_vars = xnmalloc (p, sizeof (*indep_vars));
 135   c->dependent_column = p;
 136   for (i = 0; i < p; i++)
 137     {
 138       c->indep_vars[i] = indep_vars[i];
 139     }
 140   c->indep_means = gsl_vector_alloc (p);
 141   c->indep_std = gsl_vector_alloc (p);
 142
 143   c->n_obs = n;
 144   c->n_indeps = p;
 145   c->n_coeffs = p;
 146   c->coeff = xnmalloc (p, sizeof (*c->coeff));
 147   c->cov = gsl_matrix_calloc (c->n_coeffs + 1, c->n_coeffs + 1);
 148   c->dft = n;
 149   if (!origin)
 150     c->dft--;
 151
 152   c->dfm = p;
 153   c->dfe = c->dft - c->dfm;
 154   c->intercept = 0.0;
 155   c->depvar_mean = 0.0;
 156   /*
 157      Default settings.
 158    */
 159
 160   c->refcnt = 1;
 161
 162   c->origin = origin;
 163
 164   return c;
 165 }
 166
 167
 168 void
 169 linreg_ref (struct linreg *c)
 170 {
 171   c->refcnt++;
 172 }
 173
 174 void
 175 linreg_unref (struct linreg *c)
 176 {
 177   if (--c->refcnt == 0)
 178     {
 179       gsl_vector_free (c->indep_means);
 180       gsl_vector_free (c->indep_std);
 181       gsl_matrix_free (c->cov);
 182       free (c->indep_vars);
 183       free (c->coeff);
 184       free (c);
 185     }
 186 }
 187
 188 static void
 189 post_sweep_computations (struct linreg *l, gsl_matrix *sw)
 190 {
 191   double m;
 192   size_t i;
 193   size_t j;
 194   int rc;
 195
 196   assert (sw != NULL);
 197   assert (l != NULL);
 198
 199   l->sse = gsl_matrix_get (sw, l->n_indeps, l->n_indeps);
 200   l->mse = l->sse / l->dfe;
 201   /*
 202     Get the intercept.
 203   */
 204   m = l->depvar_mean;
 205   for (i = 0; i < l->n_indeps; i++)
 206     {
 207       double tmp = gsl_matrix_get (sw, i, l->n_indeps);
 208       l->coeff[i] = tmp;
 209       m -= tmp * linreg_get_indep_variable_mean (l, i);
 210     }
 211   /*
 212     Get the covariance matrix of the parameter estimates.
 213     Only the upper triangle is necessary.
 214   */
 215
 216   /*
 217     The loops below do not compute the entries related
 218     to the estimated intercept.
 219   */
 220   for (i = 0; i < l->n_indeps; i++)
 221     for (j = i; j < l->n_indeps; j++)
 222       {
 223         double tmp = -1.0 * l->mse * gsl_matrix_get (sw, i, j);
 224         gsl_matrix_set (l->cov, i + 1, j + 1, tmp);
 225       }
 226
 227   if (! l->origin)
 228     {
 229       gsl_matrix *xm;
 230       gsl_matrix_view xtx;
 231       gsl_matrix_view xmxtx;
 232       /*
 233         Get the covariances related to the intercept.
 234       */
 235       xtx = gsl_matrix_submatrix (sw, 0, 0, l->n_indeps, l->n_indeps);
 236       xmxtx = gsl_matrix_submatrix (l->cov, 0, 1, 1, l->n_indeps);
 237       xm = gsl_matrix_calloc (1, l->n_indeps);
 238       for (i = 0; i < xm->size2; i++)
 239         {
 240           gsl_matrix_set (xm, 0, i,
 241                           linreg_get_indep_variable_mean (l, i));
 242         }
 243       rc = gsl_blas_dsymm (CblasRight, CblasUpper, l->mse,
 244                            &xtx.matrix, xm, 0.0, &xmxtx.matrix);
 245       gsl_matrix_free (xm);
 246       if (rc == GSL_SUCCESS)
 247         {
 248           double tmp = l->mse / l->n_obs;
 249           for (i = 1; i < 1 + l->n_indeps; i++)
 250             {
 251               tmp -= gsl_matrix_get (l->cov, 0, i)
 252                 * linreg_get_indep_variable_mean (l, i - 1);
 253             }
 254           gsl_matrix_set (l->cov, 0, 0, tmp);
 255
 256           l->intercept = m;
 257         }
 258       else
 259         {
 260           fprintf (stderr, "%s:%d:gsl_blas_dsymm: %s\n",
 261                    __FILE__, __LINE__, gsl_strerror (rc));
 262           exit (rc);
 263         }
 264     }
 265 }
 266
 267 /*
 268   Predict the value of the dependent variable with the new set of
 269   predictors. VALS are assumed to be in the order corresponding to the
 270   order of the coefficients in the linreg struct.
 271  */
 272 double
 273 linreg_predict (const struct linreg *c, const double *vals, size_t n_vals)
 274 {
 275   size_t j;
 276   double result;
 277
 278   assert (n_vals == c->n_coeffs);
 279   if (vals == NULL || c == NULL)
 280     {
 281       return GSL_NAN;
 282     }
 283   if (c->coeff == NULL)
 284     {
 285       /* The stupid model: just guess the mean. */
 286       return c->depvar_mean;
 287     }
 288   result = c->intercept;
 289
 290   for (j = 0; j < n_vals; j++)
 291     {
 292       result += linreg_coeff (c, j) * vals[j];
 293     }
 294
 295   return result;
 296 }
 297
 298 double
 299 linreg_residual (const struct linreg *c, double obs, const double *vals, size_t n_vals)
 300 {
 301   if (vals == NULL || c == NULL)
 302     {
 303       return GSL_NAN;
 304     }
 305   return (obs - linreg_predict (c, vals, n_vals));
 306 }
 307
 308 /*
 309   Mean of the independent variable.
 310  */
 311 double
 312 linreg_get_indep_variable_mean (const struct linreg *c, size_t j)
 313 {
 314   assert (c != NULL);
 315   return gsl_vector_get (c->indep_means, j);
 316 }
 317
 318 void
 319 linreg_set_indep_variable_mean (struct linreg *c, size_t j, double m)
 320 {
 321   assert (c != NULL);
 322   gsl_vector_set (c->indep_means, j, m);
 323 }
 324
 325 static void
 326 linreg_fit_qr (const gsl_matrix *cov, struct linreg *l)
 327 {
 328   double intcpt_coef = 0.0;
 329   double intercept_variance = 0.0;
 330   gsl_matrix *xtx;
 331   gsl_matrix *q;
 332   gsl_matrix *r;
 333   gsl_vector *xty;
 334   gsl_vector *tau;
 335   gsl_vector *params;
 336   size_t i;
 337   size_t j;
 338
 339   xtx = gsl_matrix_alloc (cov->size1 - 1, cov->size2 - 1);
 340   xty = gsl_vector_alloc (cov->size1 - 1);
 341   tau = gsl_vector_alloc (cov->size1 - 1);
 342   params = gsl_vector_alloc (cov->size1 - 1);
 343
 344   for (i = 0; i < xtx->size1; i++)
 345     {
 346       gsl_vector_set (xty, i, gsl_matrix_get (cov, cov->size2 - 1, i));
 347       for (j = 0; j < xtx->size2; j++)
 348         {
 349           gsl_matrix_set (xtx, i, j, gsl_matrix_get (cov, i, j));
 350         }
 351     }
 352   gsl_linalg_QR_decomp (xtx, tau);
 353   q = gsl_matrix_alloc (xtx->size1, xtx->size2);
 354   r = gsl_matrix_alloc (xtx->size1, xtx->size2);
 355
 356   gsl_linalg_QR_unpack (xtx, tau, q, r);
 357   gsl_linalg_QR_solve (xtx, tau, xty, params);
 358   for (i = 0; i < params->size; i++)
 359     {
 360       l->coeff[i] = gsl_vector_get (params, i);
 361     }
 362   l->sst = gsl_matrix_get (cov, cov->size1 - 1, cov->size2 - 1);
 363   l->ssm = 0.0;
 364   for (i = 0; i < l->n_indeps; i++)
 365     {
 366       l->ssm += gsl_vector_get (xty, i) * l->coeff[i];
 367     }
 368   l->sse = l->sst - l->ssm;
 369
 370   gsl_blas_dtrsm (CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit, linreg_mse (l),
 371                   r, q);
 372   /* Copy the lower triangle into the upper triangle. */
 373   for (i = 0; i < q->size1; i++)
 374     {
 375       gsl_matrix_set (l->cov, i + 1, i + 1, gsl_matrix_get (q, i, i));
 376       for (j = i + 1; j < q->size2; j++)
 377         {
 378           intercept_variance -= 2.0 * gsl_matrix_get (q, i, j) *
 379             linreg_get_indep_variable_mean (l, i) *
 380             linreg_get_indep_variable_mean (l, j);
 381           gsl_matrix_set (q, i, j, gsl_matrix_get (q, j, i));
 382         }
 383     }
 384
 385   if (!l->origin)
 386     {
 387       l->intercept = linreg_get_depvar_mean (l);
 388       for (i = 0; i < l->n_indeps; i++)
 389         {
 390           double tmp = linreg_get_indep_variable_mean (l, i);
 391           l->intercept -= l->coeff[i] * tmp;
 392           intercept_variance += tmp * tmp * gsl_matrix_get (q, i, i);
 393         }
 394
 395       /* Covariances related to the intercept. */
 396       intercept_variance += linreg_mse (l) / linreg_n_obs (l);
 397       gsl_matrix_set (l->cov, 0, 0, intercept_variance);
 398       for (i = 0; i < q->size1; i++)
 399         {
 400           for (j = 0; j < q->size2; j++)
 401             {
 402               intcpt_coef -= gsl_matrix_get (q, i, j)
 403                 * linreg_get_indep_variable_mean (l, j);
 404             }
 405           gsl_matrix_set (l->cov, 0, i + 1, intcpt_coef);
 406           gsl_matrix_set (l->cov, i + 1, 0, intcpt_coef);
 407           intcpt_coef = 0.0;
 408         }
 409     }
 410
 411   gsl_matrix_free (q);
 412   gsl_matrix_free (r);
 413   gsl_vector_free (xty);
 414   gsl_vector_free (tau);
 415   gsl_matrix_free (xtx);
 416   gsl_vector_free (params);
 417 }
 418
 419 #define REG_LARGE_DATA 1000
 420
 421 /*
 422   Estimate the model parameters from the covariance matrix. This
 423   function assumes the covariance entries corresponding to the
 424   dependent variable are in the final row and column of the covariance
 425   matrix.
 426 */
 427 void
 428 linreg_fit (const gsl_matrix *cov, struct linreg *l)
 429 {
 430   assert (l != NULL);
 431   assert (cov != NULL);
 432
 433   l->sst = gsl_matrix_get (cov, cov->size1 - 1, cov->size2 - 1);
 434
 435 #if 0
 436   /*  This QR decomposition path seems to produce the incorrect
 437       values.  See https://savannah.gnu.org/bugs/?51373  */
 438   if ((l->n_obs * l->n_obs > l->n_indeps) && (l->n_obs > REG_LARGE_DATA))
 439     {
 440       /*
 441         For large data sets, use QR decomposition.
 442       */
 443       linreg_fit_qr (cov, l);
 444     }
 445   else
 446 #endif
 447     {
 448       gsl_matrix *params = gsl_matrix_calloc (cov->size1, cov->size2);
 449       gsl_matrix_memcpy (params, cov);
 450       reg_sweep (params, l->dependent_column);
 451       post_sweep_computations (l, params);
 452       gsl_matrix_free (params);
 453     }
 454 }
 455
 456 double
 457 linreg_mse (const struct linreg *c)
 458 {
 459   assert (c != NULL);
 460   return (c->sse / c->dfe);
 461 }
 462
 463 double
 464 linreg_intercept (const struct linreg *c)
 465 {
 466   return c->intercept;
 467 }
 468
 469 const gsl_matrix *
 470 linreg_cov (const struct linreg *c)
 471 {
 472   return c->cov;
 473 }
 474
 475 double
 476 linreg_coeff (const struct linreg *c, size_t i)
 477 {
 478   return (c->coeff[i]);
 479 }
 480
 481 const struct variable *
 482 linreg_indep_var (const struct linreg *c, size_t i)
 483 {
 484   return (c->indep_vars[i]);
 485 }
 486
 487 int
 488 linreg_n_indeps (const struct linreg *c)
 489 {
 490   return c->n_indeps;
 491 }
 492
 493
 494 const struct variable *
 495 linreg_dep_var (const struct linreg *c)
 496 {
 497   return c->depvar;
 498 }
 499
 500
 501 size_t
 502 linreg_n_coeffs (const struct linreg *c)
 503 {
 504   return c->n_coeffs;
 505 }
 506
 507 double
 508 linreg_n_obs (const struct linreg *c)
 509 {
 510   return c->n_obs;
 511 }
 512
 513 double
 514 linreg_sse (const struct linreg *c)
 515 {
 516   return c->sse;
 517 }
 518
 519 double
 520 linreg_ssreg (const struct linreg *c)
 521 {
 522   return (c->sst - c->sse);
 523 }
 524
 525 double linreg_sst (const struct linreg *c)
 526 {
 527   return c->sst;
 528 }
 529
 530 double
 531 linreg_dfmodel ( const struct linreg *c)
 532 {
 533   return c->dfm;
 534 }
 535
 536 double
 537 linreg_dferror ( const struct linreg *c)
 538 {
 539   return c->dfe;
 540 }
 541
 542 double
 543 linreg_dftotal ( const struct linreg *c)
 544 {
 545   return c->dft;
 546 }
 547
 548 void
 549 linreg_set_depvar_mean (struct linreg *c, double x)
 550 {
 551   c->depvar_mean = x;
 552 }
 553
 554 double
 555 linreg_get_depvar_mean (const struct linreg *c)
 556 {
 557   return c->depvar_mean;
 558 }