src/math/linreg.c

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