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

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2005, 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/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
  61 const struct variable **
  62 linreg_get_vars (const linreg *c)
  63 {
  64   return c->indep_vars;
  65 }
  66
  67 /*
  68   Allocate a linreg and return a pointer to it. n is the number of
  69   cases, p is the number of independent variables.
  70  */
  71 linreg *
  72 linreg_alloc (const struct variable *depvar, const struct variable **indep_vars,
  73               double n, size_t p)
  74 {
  75   linreg *c;
  76   size_t i;
  77
  78   c = xmalloc (sizeof (*c));
  79   c->depvar = depvar;
  80   c->indep_vars = xnmalloc (p, sizeof (*indep_vars));
  81   c->dependent_column = p;
  82   for (i = 0; i < p; i++)
  83     {
  84       c->indep_vars[i] = indep_vars[i];
  85     }
  86   c->indep_means = gsl_vector_alloc (p);
  87   c->indep_std = gsl_vector_alloc (p);
  88
  89   c->n_obs = n;
  90   c->n_indeps = p;
  91   c->n_coeffs = p;
  92   c->coeff = xnmalloc (p, sizeof (*c->coeff));
  93   c->cov = gsl_matrix_calloc (c->n_coeffs + 1, c->n_coeffs + 1);
  94   c->dft = n - 1;
  95   c->dfm = p;
  96   c->dfe = c->dft - c->dfm;
  97   c->intercept = 0.0;
  98   c->depvar_mean = 0.0;
  99   /*
 100      Default settings.
 101    */
 102   c->method = LINREG_SWEEP;
 103
 104   c->refcnt = 1;
 105
 106   return c;
 107 }
 108
 109
 110 void
 111 linreg_ref (linreg *c)
 112 {
 113   c->refcnt++;
 114 }
 115
 116 void
 117 linreg_unref (linreg *c)
 118 {
 119   if (--c->refcnt == 0)
 120     {
 121       gsl_vector_free (c->indep_means);
 122       gsl_vector_free (c->indep_std);
 123       gsl_matrix_free (c->cov);
 124       free (c->indep_vars);
 125       free (c->coeff);
 126       free (c);
 127     }
 128 }
 129
 130 static void
 131 post_sweep_computations (linreg *l, gsl_matrix *sw)
 132 {
 133   gsl_matrix *xm;
 134   gsl_matrix_view xtx;
 135   gsl_matrix_view xmxtx;
 136   double m;
 137   double tmp;
 138   size_t i;
 139   size_t j;
 140   int rc;
 141
 142   assert (sw != NULL);
 143   assert (l != NULL);
 144
 145   l->sse = gsl_matrix_get (sw, l->n_indeps, l->n_indeps);
 146   l->mse = l->sse / l->dfe;
 147   /*
 148     Get the intercept.
 149   */
 150   m = l->depvar_mean;
 151   for (i = 0; i < l->n_indeps; i++)
 152     {
 153       tmp = gsl_matrix_get (sw, i, l->n_indeps);
 154       l->coeff[i] = tmp;
 155       m -= tmp * linreg_get_indep_variable_mean (l, i);
 156     }
 157   /*
 158     Get the covariance matrix of the parameter estimates.
 159     Only the upper triangle is necessary.
 160   */
 161
 162   /*
 163     The loops below do not compute the entries related
 164     to the estimated intercept.
 165   */
 166   for (i = 0; i < l->n_indeps; i++)
 167     for (j = i; j < l->n_indeps; j++)
 168       {
 169         tmp = -1.0 * l->mse * gsl_matrix_get (sw, i, j);
 170         gsl_matrix_set (l->cov, i + 1, j + 1, tmp);
 171       }
 172   /*
 173     Get the covariances related to the intercept.
 174   */
 175   xtx = gsl_matrix_submatrix (sw, 0, 0, l->n_indeps, l->n_indeps);
 176   xmxtx = gsl_matrix_submatrix (l->cov, 0, 1, 1, l->n_indeps);
 177   xm = gsl_matrix_calloc (1, l->n_indeps);
 178   for (i = 0; i < xm->size2; i++)
 179     {
 180       gsl_matrix_set (xm, 0, i,
 181                       linreg_get_indep_variable_mean (l, i));
 182     }
 183   rc = gsl_blas_dsymm (CblasRight, CblasUpper, l->mse,
 184                        &xtx.matrix, xm, 0.0, &xmxtx.matrix);
 185   gsl_matrix_free (xm);
 186   if (rc == GSL_SUCCESS)
 187     {
 188       tmp = l->mse / l->n_obs;
 189       for (i = 1; i < 1 + l->n_indeps; i++)
 190         {
 191           tmp -= gsl_matrix_get (l->cov, 0, i)
 192             * linreg_get_indep_variable_mean (l, i - 1);
 193         }
 194       gsl_matrix_set (l->cov, 0, 0, tmp);
 195
 196       l->intercept = m;
 197     }
 198   else
 199     {
 200       fprintf (stderr, "%s:%d:gsl_blas_dsymm: %s\n",
 201                __FILE__, __LINE__, gsl_strerror (rc));
 202       exit (rc);
 203     }
 204 }
 205
 206 /*
 207   Predict the value of the dependent variable with the new set of
 208   predictors. VALS are assumed to be in the order corresponding to the
 209   order of the coefficients in the linreg struct.
 210  */
 211 double
 212 linreg_predict (const linreg *c, const double *vals, size_t n_vals)
 213 {
 214   size_t j;
 215   double result;
 216
 217   assert (n_vals = c->n_coeffs);
 218   if (vals == NULL || c == NULL)
 219     {
 220       return GSL_NAN;
 221     }
 222   if (c->coeff == NULL)
 223     {
 224       /* The stupid model: just guess the mean. */
 225       return c->depvar_mean;
 226     }
 227   result = c->intercept;
 228
 229   for (j = 0; j < n_vals; j++)
 230     {
 231       result += linreg_coeff (c, j) * vals[j];
 232     }
 233
 234   return result;
 235 }
 236
 237 double
 238 linreg_residual (const linreg *c, double obs, const double *vals, size_t n_vals)
 239 {
 240   if (vals == NULL || c == NULL)
 241     {
 242       return GSL_NAN;
 243     }
 244   return (obs - linreg_predict (c, vals, n_vals));
 245 }
 246
 247 /*
 248   Mean of the independent variable.
 249  */
 250 double linreg_get_indep_variable_mean (const linreg *c, size_t j)
 251 {
 252   assert (c != NULL);
 253   return gsl_vector_get (c->indep_means, j);
 254 }
 255
 256 void linreg_set_indep_variable_mean (linreg *c, size_t j, double m)
 257 {
 258   assert (c != NULL);
 259   gsl_vector_set (c->indep_means, j, m);
 260 }
 261
 262 static void
 263 linreg_fit_qr (const gsl_matrix *cov, linreg *l)
 264 {
 265   double intcpt_coef = 0.0;
 266   double intercept_variance = 0.0;
 267   gsl_matrix *xtx;
 268   gsl_matrix *q;
 269   gsl_matrix *r;
 270   gsl_vector *xty;
 271   gsl_vector *tau;
 272   gsl_vector *params;
 273   double tmp = 0.0;
 274   size_t i;
 275   size_t j;
 276
 277   xtx = gsl_matrix_alloc (cov->size1 - 1, cov->size2 - 1);
 278   xty = gsl_vector_alloc (cov->size1 - 1);
 279   tau = gsl_vector_alloc (cov->size1 - 1);
 280   params = gsl_vector_alloc (cov->size1 - 1);
 281
 282   for (i = 0; i < xtx->size1; i++)
 283     {
 284       gsl_vector_set (xty, i, gsl_matrix_get (cov, cov->size2 - 1, i));
 285       for (j = 0; j < xtx->size2; j++)
 286         {
 287           gsl_matrix_set (xtx, i, j, gsl_matrix_get (cov, i, j));
 288         }
 289     }
 290   gsl_linalg_QR_decomp (xtx, tau);
 291   q = gsl_matrix_alloc (xtx->size1, xtx->size2);
 292   r = gsl_matrix_alloc (xtx->size1, xtx->size2);
 293
 294   gsl_linalg_QR_unpack (xtx, tau, q, r);
 295   gsl_linalg_QR_solve (xtx, tau, xty, params);
 296   for (i = 0; i < params->size; i++)
 297     {
 298       l->coeff[i] = gsl_vector_get (params, i);
 299     }
 300   l->sst = gsl_matrix_get (cov, cov->size1 - 1, cov->size2 - 1);
 301   l->ssm = 0.0;
 302   for (i = 0; i < l->n_indeps; i++)
 303     {
 304       l->ssm += gsl_vector_get (xty, i) * l->coeff[i];
 305     }
 306   l->sse = l->sst - l->ssm;
 307
 308   gsl_blas_dtrsm (CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit, linreg_mse (l),
 309                   r, q);
 310   /* Copy the lower triangle into the upper triangle. */
 311   for (i = 0; i < q->size1; i++)
 312     {
 313       gsl_matrix_set (l->cov, i + 1, i + 1, gsl_matrix_get (q, i, i));
 314       for (j = i + 1; j < q->size2; j++)
 315         {
 316           intercept_variance -= 2.0 * gsl_matrix_get (q, i, j) *
 317             linreg_get_indep_variable_mean (l, i) *
 318             linreg_get_indep_variable_mean (l, j);
 319           gsl_matrix_set (q, i, j, gsl_matrix_get (q, j, i));
 320         }
 321     }
 322   l->intercept = linreg_get_depvar_mean (l);
 323   tmp = 0.0;
 324   for (i = 0; i < l->n_indeps; i++)
 325     {
 326       tmp = linreg_get_indep_variable_mean (l, i);
 327       l->intercept -= l->coeff[i] * tmp;
 328       intercept_variance += tmp * tmp * gsl_matrix_get (q, i, i);
 329     }
 330
 331   /* Covariances related to the intercept. */
 332   intercept_variance += linreg_mse (l) / linreg_n_obs (l);
 333   gsl_matrix_set (l->cov, 0, 0, intercept_variance);
 334   for (i = 0; i < q->size1; i++)
 335     {
 336       for (j = 0; j < q->size2; j++)
 337         {
 338           intcpt_coef -= gsl_matrix_get (q, i, j)
 339             * linreg_get_indep_variable_mean (l, j);
 340         }
 341       gsl_matrix_set (l->cov, 0, i + 1, intcpt_coef);
 342       gsl_matrix_set (l->cov, i + 1, 0, intcpt_coef);
 343       intcpt_coef = 0.0;
 344     }
 345
 346   gsl_matrix_free (q);
 347   gsl_matrix_free (r);
 348   gsl_vector_free (xty);
 349   gsl_vector_free (tau);
 350   gsl_matrix_free (xtx);
 351   gsl_vector_free (params);
 352 }
 353
 354 /*
 355   Estimate the model parameters from the covariance matrix. This
 356   function assumes the covariance entries corresponding to the
 357   dependent variable are in the final row and column of the covariance
 358   matrix.
 359 */
 360 void
 361 linreg_fit (const gsl_matrix *cov, linreg *l)
 362 {
 363   assert (l != NULL);
 364   assert (cov != NULL);
 365
 366   l->sst = gsl_matrix_get (cov, cov->size1 - 1, cov->size2 - 1);
 367   if (l->method == LINREG_SWEEP)
 368     {
 369       gsl_matrix *params;
 370       params = gsl_matrix_calloc (cov->size1, cov->size2);
 371       gsl_matrix_memcpy (params, cov);
 372       reg_sweep (params, l->dependent_column);
 373       post_sweep_computations (l, params);
 374       gsl_matrix_free (params);
 375     }
 376   else if (l->method == LINREG_QR)
 377     {
 378       linreg_fit_qr (cov, l);
 379     }
 380 }
 381
 382 double linreg_mse (const linreg *c)
 383 {
 384   assert (c != NULL);
 385   return (c->sse / c->dfe);
 386 }
 387
 388 double linreg_intercept (const linreg *c)
 389 {
 390   return c->intercept;
 391 }
 392
 393 const gsl_matrix *
 394 linreg_cov (const linreg *c)
 395 {
 396   return c->cov;
 397 }
 398
 399 double
 400 linreg_coeff (const linreg *c, size_t i)
 401 {
 402   return (c->coeff[i]);
 403 }
 404
 405 const struct variable *
 406 linreg_indep_var (const linreg *c, size_t i)
 407 {
 408   return (c->indep_vars[i]);
 409 }
 410
 411 size_t
 412 linreg_n_coeffs (const linreg *c)
 413 {
 414   return c->n_coeffs;
 415 }
 416
 417 double
 418 linreg_n_obs (const linreg *c)
 419 {
 420   return c->n_obs;
 421 }
 422
 423 double
 424 linreg_sse (const linreg *c)
 425 {
 426   return c->sse;
 427 }
 428
 429 double
 430 linreg_ssreg (const linreg *c)
 431 {
 432   return (c->sst - c->sse);
 433 }
 434
 435 double linreg_sst (const linreg *c)
 436 {
 437   return c->sst;
 438 }
 439
 440 double
 441 linreg_dfmodel ( const linreg *c)
 442 {
 443   return c->dfm;
 444 }
 445
 446 void
 447 linreg_set_depvar_mean (linreg *c, double x)
 448 {
 449   c->depvar_mean = x;
 450 }
 451
 452 double
 453 linreg_get_depvar_mean (const linreg *c)
 454 {
 455   return c->depvar_mean;
 456 }