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   c->ssx = gsl_vector_alloc (p);        /* Sums of squares for the
  85                                            independent variables.
  86                                          */
  87   c->ss_indeps = gsl_vector_alloc (p);  /* Sums of squares for the
  88                                            model parameters.
  89                                          */
  90   c->n_obs = n;
  91   c->n_indeps = p;
  92   c->n_coeffs = p;
  93   c->coeff = xnmalloc (p, sizeof (*c->coeff));
  94   c->cov = gsl_matrix_alloc (c->n_coeffs + 1, c->n_coeffs + 1);
  95   c->dft = n - 1;
  96   c->dfm = p;
  97   c->dfe = c->dft - c->dfm;
  98   c->intercept = 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       gsl_vector_free (c->ssx);
 119       free (c->indep_vars);
 120       free (c->coeff);
 121       free (c);
 122     }
 123   return true;
 124 }
 125
 126 static void
 127 post_sweep_computations (linreg *l, gsl_matrix *sw)
 128 {
 129   gsl_matrix *xm;
 130   gsl_matrix_view xtx;
 131   gsl_matrix_view xmxtx;
 132   double m;
 133   double tmp;
 134   size_t i;
 135   size_t j;
 136   int rc;
 137
 138   assert (sw != NULL);
 139   assert (l != NULL);
 140
 141   l->sse = gsl_matrix_get (sw, l->n_indeps, l->n_indeps);
 142   l->mse = l->sse / l->dfe;
 143   /*
 144     Get the intercept.
 145   */
 146   m = l->depvar_mean;
 147   for (i = 0; i < l->n_indeps; i++)
 148     {
 149       tmp = gsl_matrix_get (sw, i, l->n_indeps);
 150       l->coeff[i] = tmp;
 151       m -= tmp * linreg_get_indep_variable_mean (l, i);
 152     }
 153   /*
 154     Get the covariance matrix of the parameter estimates.
 155     Only the upper triangle is necessary.
 156   */
 157
 158   /*
 159     The loops below do not compute the entries related
 160     to the estimated intercept.
 161   */
 162   for (i = 0; i < l->n_indeps; i++)
 163     for (j = i; j < l->n_indeps; j++)
 164       {
 165         tmp = -1.0 * l->mse * gsl_matrix_get (sw, i, j);
 166         gsl_matrix_set (l->cov, i + 1, j + 1, tmp);
 167       }
 168   /*
 169     Get the covariances related to the intercept.
 170   */
 171   xtx = gsl_matrix_submatrix (sw, 0, 0, l->n_indeps, l->n_indeps);
 172   xmxtx = gsl_matrix_submatrix (l->cov, 0, 1, 1, l->n_indeps);
 173   xm = gsl_matrix_calloc (1, l->n_indeps);
 174   for (i = 0; i < xm->size2; i++)
 175     {
 176       gsl_matrix_set (xm, 0, i,
 177                       linreg_get_indep_variable_mean (l, i));
 178     }
 179   rc = gsl_blas_dsymm (CblasRight, CblasUpper, l->mse,
 180                        &xtx.matrix, xm, 0.0, &xmxtx.matrix);
 181   gsl_matrix_free (xm);
 182   if (rc == GSL_SUCCESS)
 183     {
 184       tmp = l->mse / l->n_obs;
 185       for (i = 1; i < 1 + l->n_indeps; i++)
 186         {
 187           tmp -= gsl_matrix_get (l->cov, 0, i)
 188             * linreg_get_indep_variable_mean (l, i - 1);
 189         }
 190       gsl_matrix_set (l->cov, 0, 0, tmp);
 191
 192       l->intercept = m;
 193     }
 194   else
 195     {
 196       fprintf (stderr, "%s:%d:gsl_blas_dsymm: %s\n",
 197                __FILE__, __LINE__, gsl_strerror (rc));
 198       exit (rc);
 199     }
 200 }
 201
 202 /*
 203   Predict the value of the dependent variable with the new set of
 204   predictors. VALS are assumed to be in the order corresponding to the
 205   order of the coefficients in the linreg struct.
 206  */
 207 double
 208 linreg_predict (const linreg *c, const double *vals, size_t n_vals)
 209 {
 210   size_t j;
 211   double result;
 212
 213   assert (n_vals = c->n_coeffs);
 214   if (vals == NULL || c == NULL)
 215     {
 216       return GSL_NAN;
 217     }
 218   if (c->coeff == NULL)
 219     {
 220       /* The stupid model: just guess the mean. */
 221       return c->depvar_mean;
 222     }
 223   result = c->intercept;
 224
 225   for (j = 0; j < n_vals; j++)
 226     {
 227       result += linreg_coeff (c, j) * vals[j];
 228     }
 229
 230   return result;
 231 }
 232
 233 double
 234 linreg_residual (const linreg *c, double obs, const double *vals, size_t n_vals)
 235 {
 236   if (vals == NULL || c == NULL)
 237     {
 238       return GSL_NAN;
 239     }
 240   return (obs - linreg_predict (c, vals, n_vals));
 241 }
 242
 243 double linreg_get_indep_variable_sd (linreg *c, size_t j)
 244 {
 245   assert (c != NULL);
 246   return gsl_vector_get (c->indep_std, j);
 247 }
 248
 249 void linreg_set_indep_variable_sd (linreg *c, size_t j, double s)
 250 {
 251   assert (c != NULL);
 252   gsl_vector_set (c->indep_std, j, s);
 253 }
 254
 255 /*
 256   Mean of the independent variable.
 257  */
 258 double linreg_get_indep_variable_mean (linreg *c, size_t j)
 259 {
 260   assert (c != NULL);
 261   return gsl_vector_get (c->indep_means, j);
 262 }
 263
 264 void linreg_set_indep_variable_mean (linreg *c, size_t j, double m)
 265 {
 266   assert (c != NULL);
 267   gsl_vector_set (c->indep_means, j, m);
 268 }
 269
 270 static void
 271 linreg_fit_qr (const gsl_matrix *cov, linreg *l)
 272 {
 273   gsl_matrix *xtx;
 274   gsl_matrix *q;
 275   gsl_matrix *r;
 276   gsl_vector *xty;
 277   gsl_vector *tau;
 278   gsl_vector *params;
 279   double tmp = 0.0;
 280   size_t i;
 281   size_t j;
 282
 283   xtx = gsl_matrix_alloc (cov->size1 - 1, cov->size2 - 1);
 284   xty = gsl_vector_alloc (cov->size1 - 1);
 285   tau = gsl_vector_alloc (cov->size1 - 1);
 286   params = gsl_vector_alloc (cov->size1 - 1);
 287
 288   for (i = 0; i < xtx->size1; i++)
 289     {
 290       gsl_vector_set (xty, i, gsl_matrix_get (cov, cov->size2 - 1, i));
 291       for (j = 0; j < xtx->size2; j++)
 292         {
 293           gsl_matrix_set (xtx, i, j, gsl_matrix_get (cov, i, j));
 294         }
 295     }
 296   gsl_linalg_QR_decomp (xtx, tau);
 297   q = gsl_matrix_alloc (xtx->size1, xtx->size2);
 298   r = gsl_matrix_alloc (xtx->size1, xtx->size2);
 299   gsl_linalg_QR_unpack (xtx, tau, q, r);
 300   gsl_linalg_QR_solve (xtx, tau, xty, params);
 301   for (i = 0; i < params->size; i++)
 302     {
 303       l->coeff[i] = gsl_vector_get (params, i);
 304     }
 305
 306   l->intercept = l->depvar_mean;
 307   for (i = 0; i < l->n_indeps; i++)
 308     {
 309       l->intercept -= l->coeff[i] * linreg_get_indep_variable_mean (l, i);
 310     }
 311
 312   l->sse = gsl_matrix_get (cov, cov->size1 - 1, cov->size2 - 1);
 313   for (i = 0; i < l->n_indeps; i++)
 314     {
 315       tmp += gsl_vector_get (xty, i) * l->coeff[i];
 316     }
 317   l->sse -= 2.0 * tmp;
 318   for (i = 0; i < xtx->size1; i++)
 319     {
 320       tmp = 0.0;
 321       for (j = i; j < xtx->size2; j++)
 322         {
 323           tmp += gsl_matrix_get (xtx, i, j) * l->coeff[j];
 324         }
 325       l->sse += tmp * tmp;
 326     }
 327
 328 #if 0
 329   p = l->hat->size1 - 1;
 330   for (i = 0; i < l->cov->size1 - 1; i++)
 331     {
 332       gsl_matrix_set (l->cov, p - i, p - i, 1.0 / gsl_matrix_get (r, p - i, p - i));
 333       for (j = 0; j < i; j++)
 334         {
 335           tmp = -1.0 * gsl_matrix_get (r, p - i, p - j);
 336           tmp /= gsl_matrix_get (r, p - i, p - i) * gsl_matrix_get (r, p - j, p - j);
 337           gsl_matrix_set (l->cov, p - i, p - j, tmp);
 338         }
 339     }
 340 #endif
 341   gsl_matrix_transpose_memcpy (l->cov, q);
 342   gsl_blas_dtrsm (CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, 1.0,
 343                   r, l->cov);
 344
 345   gsl_matrix_free (q);
 346   gsl_matrix_free (r);
 347   gsl_vector_free (xty);
 348   gsl_vector_free (tau);
 349   gsl_matrix_free (xtx);
 350   gsl_vector_free (params);
 351 }
 352
 353 /*
 354   Estimate the model parameters from the covariance matrix. This
 355   function assumes the covariance entries corresponding to the
 356   dependent variable are in the final row and column of the covariance
 357   matrix.
 358 */
 359 void
 360 linreg_fit (const gsl_matrix *cov, linreg *l)
 361 {
 362   gsl_matrix *params;
 363   assert (l != NULL);
 364   assert (cov != NULL);
 365
 366   params = gsl_matrix_calloc (cov->size1, cov->size2);
 367   gsl_matrix_memcpy (params, cov);
 368   l->sst = gsl_matrix_get (cov, cov->size1 - 1, cov->size2 - 1);
 369
 370   if (l->method == LINREG_SWEEP)
 371     {
 372       reg_sweep (params);
 373       post_sweep_computations (l, params);
 374     }
 375   else if (l->method == LINREG_QR)
 376     {
 377       linreg_fit_qr (params, l);
 378     }
 379   gsl_matrix_free (params);
 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 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 size_t
 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 }