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 linreg_get_indep_variable_mean (const struct linreg *c, size_t j)
 312 {
 313   assert (c != NULL);
 314   return gsl_vector_get (c->indep_means, j);
 315 }
 316
 317 void linreg_set_indep_variable_mean (struct linreg *c, size_t j, double m)
 318 {
 319   assert (c != NULL);
 320   gsl_vector_set (c->indep_means, j, m);
 321 }
 322
 323 static void
 324 linreg_fit_qr (const gsl_matrix *cov, struct linreg *l)
 325 {
 326   double intcpt_coef = 0.0;
 327   double intercept_variance = 0.0;
 328   gsl_matrix *xtx;
 329   gsl_matrix *q;
 330   gsl_matrix *r;
 331   gsl_vector *xty;
 332   gsl_vector *tau;
 333   gsl_vector *params;
 334   size_t i;
 335   size_t j;
 336
 337   xtx = gsl_matrix_alloc (cov->size1 - 1, cov->size2 - 1);
 338   xty = gsl_vector_alloc (cov->size1 - 1);
 339   tau = gsl_vector_alloc (cov->size1 - 1);
 340   params = gsl_vector_alloc (cov->size1 - 1);
 341
 342   for (i = 0; i < xtx->size1; i++)
 343     {
 344       gsl_vector_set (xty, i, gsl_matrix_get (cov, cov->size2 - 1, i));
 345       for (j = 0; j < xtx->size2; j++)
 346         {
 347           gsl_matrix_set (xtx, i, j, gsl_matrix_get (cov, i, j));
 348         }
 349     }
 350   gsl_linalg_QR_decomp (xtx, tau);
 351   q = gsl_matrix_alloc (xtx->size1, xtx->size2);
 352   r = gsl_matrix_alloc (xtx->size1, xtx->size2);
 353
 354   gsl_linalg_QR_unpack (xtx, tau, q, r);
 355   gsl_linalg_QR_solve (xtx, tau, xty, params);
 356   for (i = 0; i < params->size; i++)
 357     {
 358       l->coeff[i] = gsl_vector_get (params, i);
 359     }
 360   l->sst = gsl_matrix_get (cov, cov->size1 - 1, cov->size2 - 1);
 361   l->ssm = 0.0;
 362   for (i = 0; i < l->n_indeps; i++)
 363     {
 364       l->ssm += gsl_vector_get (xty, i) * l->coeff[i];
 365     }
 366   l->sse = l->sst - l->ssm;
 367
 368   gsl_blas_dtrsm (CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit, linreg_mse (l),
 369                   r, q);
 370   /* Copy the lower triangle into the upper triangle. */
 371   for (i = 0; i < q->size1; i++)
 372     {
 373       gsl_matrix_set (l->cov, i + 1, i + 1, gsl_matrix_get (q, i, i));
 374       for (j = i + 1; j < q->size2; j++)
 375         {
 376           intercept_variance -= 2.0 * gsl_matrix_get (q, i, j) *
 377             linreg_get_indep_variable_mean (l, i) *
 378             linreg_get_indep_variable_mean (l, j);
 379           gsl_matrix_set (q, i, j, gsl_matrix_get (q, j, i));
 380         }
 381     }
 382   l->intercept = linreg_get_depvar_mean (l);
 383   for (i = 0; i < l->n_indeps; i++)
 384     {
 385       double tmp = linreg_get_indep_variable_mean (l, i);
 386       l->intercept -= l->coeff[i] * tmp;
 387       intercept_variance += tmp * tmp * gsl_matrix_get (q, i, i);
 388     }
 389
 390   /* Covariances related to the intercept. */
 391   intercept_variance += linreg_mse (l) / linreg_n_obs (l);
 392   gsl_matrix_set (l->cov, 0, 0, intercept_variance);
 393   for (i = 0; i < q->size1; i++)
 394     {
 395       for (j = 0; j < q->size2; j++)
 396         {
 397           intcpt_coef -= gsl_matrix_get (q, i, j)
 398             * linreg_get_indep_variable_mean (l, j);
 399         }
 400       gsl_matrix_set (l->cov, 0, i + 1, intcpt_coef);
 401       gsl_matrix_set (l->cov, i + 1, 0, intcpt_coef);
 402       intcpt_coef = 0.0;
 403     }
 404
 405   gsl_matrix_free (q);
 406   gsl_matrix_free (r);
 407   gsl_vector_free (xty);
 408   gsl_vector_free (tau);
 409   gsl_matrix_free (xtx);
 410   gsl_vector_free (params);
 411 }
 412
 413 #define REG_LARGE_DATA 1000
 414
 415 /*
 416   Estimate the model parameters from the covariance matrix. This
 417   function assumes the covariance entries corresponding to the
 418   dependent variable are in the final row and column of the covariance
 419   matrix.
 420 */
 421 void
 422 linreg_fit (const gsl_matrix *cov, struct linreg *l)
 423 {
 424   assert (l != NULL);
 425   assert (cov != NULL);
 426
 427   l->sst = gsl_matrix_get (cov, cov->size1 - 1, cov->size2 - 1);
 428
 429   if ((l->n_obs * l->n_obs > l->n_indeps) && ( l->n_obs > REG_LARGE_DATA))
 430     {
 431       /*
 432         For large data sets, use QR decomposition.
 433       */
 434       linreg_fit_qr (cov, l);
 435     }
 436   else
 437     {
 438       gsl_matrix *params = gsl_matrix_calloc (cov->size1, cov->size2);
 439       gsl_matrix_memcpy (params, cov);
 440       reg_sweep (params, l->dependent_column);
 441       post_sweep_computations (l, params);
 442       gsl_matrix_free (params);
 443     }
 444 }
 445
 446 double linreg_mse (const struct linreg *c)
 447 {
 448   assert (c != NULL);
 449   return (c->sse / c->dfe);
 450 }
 451
 452 double linreg_intercept (const struct linreg *c)
 453 {
 454   return c->intercept;
 455 }
 456
 457 const gsl_matrix *
 458 linreg_cov (const struct linreg *c)
 459 {
 460   return c->cov;
 461 }
 462
 463 double
 464 linreg_coeff (const struct linreg *c, size_t i)
 465 {
 466   return (c->coeff[i]);
 467 }
 468
 469 const struct variable *
 470 linreg_indep_var (const struct linreg *c, size_t i)
 471 {
 472   return (c->indep_vars[i]);
 473 }
 474
 475 int
 476 linreg_n_indeps (const struct linreg *c)
 477 {
 478   return c->n_indeps;
 479 }
 480
 481
 482 const struct variable *
 483 linreg_dep_var (const struct linreg *c)
 484 {
 485   return c->depvar;
 486 }
 487
 488
 489 size_t
 490 linreg_n_coeffs (const struct linreg *c)
 491 {
 492   return c->n_coeffs;
 493 }
 494
 495 double
 496 linreg_n_obs (const struct linreg *c)
 497 {
 498   return c->n_obs;
 499 }
 500
 501 double
 502 linreg_sse (const struct linreg *c)
 503 {
 504   return c->sse;
 505 }
 506
 507 double
 508 linreg_ssreg (const struct linreg *c)
 509 {
 510   return (c->sst - c->sse);
 511 }
 512
 513 double linreg_sst (const struct linreg *c)
 514 {
 515   return c->sst;
 516 }
 517
 518 double
 519 linreg_dfmodel ( const struct linreg *c)
 520 {
 521   return c->dfm;
 522 }
 523
 524 double
 525 linreg_dferror ( const struct linreg *c)
 526 {
 527   return c->dfe;
 528 }
 529
 530 double
 531 linreg_dftotal ( const struct linreg *c)
 532 {
 533   return c->dft;
 534 }
 535
 536 void
 537 linreg_set_depvar_mean (struct linreg *c, double x)
 538 {
 539   c->depvar_mean = x;
 540 }
 541
 542 double
 543 linreg_get_depvar_mean (const struct linreg *c)
 544 {
 545   return c->depvar_mean;
 546 }