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