pintos-os.org Git - pspp/blob - 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
 268 static void
 269 linreg_fit_qr (const gsl_matrix *cov, linreg *l)
 270 {
 271   double intcpt_coef = 0.0;
 272   double intercept_variance = 0.0;
 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
 300   gsl_linalg_QR_unpack (xtx, tau, q, r);
 301   gsl_linalg_QR_solve (xtx, tau, xty, params);
 302   for (i = 0; i < params->size; i++)
 303     {
 304       l->coeff[i] = gsl_vector_get (params, i);
 305     }
 306   l->sst = gsl_matrix_get (cov, cov->size1 - 1, cov->size2 - 1);
 307   l->ssm = 0.0;
 308   for (i = 0; i < l->n_indeps; i++)
 309     {
 310       l->ssm += gsl_vector_get (xty, i) * l->coeff[i];
 311     }
 312   l->sse = l->sst - l->ssm;
 313
 314   gsl_blas_dtrsm (CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit, linreg_mse (l),
 315                   r, q);
 316   /* Copy the lower triangle into the upper triangle. */
 317   for (i = 0; i < q->size1; i++)
 318     {
 319       gsl_matrix_set (l->cov, i + 1, i + 1, gsl_matrix_get (q, i, i));
 320       for (j = i + 1; j < q->size2; j++)
 321         {
 322           intercept_variance -= 2.0 * gsl_matrix_get (q, i, j) *
 323             linreg_get_indep_variable_mean (l, i) *
 324             linreg_get_indep_variable_mean (l, j);
 325           gsl_matrix_set (q, i, j, gsl_matrix_get (q, j, i));
 326         }
 327     }
 328   l->intercept = linreg_get_depvar_mean (l);
 329   tmp = 0.0;
 330   for (i = 0; i < l->n_indeps; i++)
 331     {
 332       tmp = linreg_get_indep_variable_mean (l, i);
 333       l->intercept -= l->coeff[i] * tmp;
 334       intercept_variance += tmp * tmp * gsl_matrix_get (q, i, i);
 335     }
 336
 337   /* Covariances related to the intercept. */
 338   intercept_variance += linreg_mse (l) / linreg_n_obs (l);
 339   gsl_matrix_set (l->cov, 0, 0, intercept_variance);
 340   for (i = 0; i < q->size1; i++)
 341     {
 342       for (j = 0; j < q->size2; j++)
 343         {
 344           intcpt_coef -= gsl_matrix_get (q, i, j)
 345             * linreg_get_indep_variable_mean (l, j);
 346         }
 347       gsl_matrix_set (l->cov, 0, i + 1, intcpt_coef);
 348       gsl_matrix_set (l->cov, i + 1, 0, intcpt_coef);
 349       intcpt_coef = 0.0;
 350     }
 351
 352   gsl_matrix_free (q);
 353   gsl_matrix_free (r);
 354   gsl_vector_free (xty);
 355   gsl_vector_free (tau);
 356   gsl_matrix_free (xtx);
 357   gsl_vector_free (params);
 358 }
 359
 360 /*
 361   Estimate the model parameters from the covariance matrix. This
 362   function assumes the covariance entries corresponding to the
 363   dependent variable are in the final row and column of the covariance
 364   matrix.
 365 */
 366 void
 367 linreg_fit (const gsl_matrix *cov, linreg *l)
 368 {
 369   assert (l != NULL);
 370   assert (cov != NULL);
 371
 372   l->sst = gsl_matrix_get (cov, cov->size1 - 1, cov->size2 - 1);
 373   if (l->method == LINREG_SWEEP)
 374     {
 375       gsl_matrix *params;
 376       params = gsl_matrix_calloc (cov->size1, cov->size2);
 377       gsl_matrix_memcpy (params, cov);
 378       reg_sweep (params);
 379       post_sweep_computations (l, params);
 380       gsl_matrix_free (params);
 381     }
 382   else if (l->method == LINREG_QR)
 383     {
 384       linreg_fit_qr (cov, l);
 385     }
 386 }
 387
 388 double linreg_mse (const linreg *c)
 389 {
 390   assert (c != NULL);
 391   return (c->sse / c->dfe);
 392 }
 393
 394 double linreg_intercept (const linreg *c)
 395 {
 396   return c->intercept;
 397 }
 398
 399 gsl_matrix *
 400 linreg_cov (const linreg *c)
 401 {
 402   return c->cov;
 403 }
 404
 405 double
 406 linreg_coeff (const linreg *c, size_t i)
 407 {
 408   return (c->coeff[i]);
 409 }
 410
 411 const struct variable *
 412 linreg_indep_var (const linreg *c, size_t i)
 413 {
 414   return (c->indep_vars[i]);
 415 }
 416
 417 size_t
 418 linreg_n_coeffs (const linreg *c)
 419 {
 420   return c->n_coeffs;
 421 }
 422
 423 double
 424 linreg_n_obs (const linreg *c)
 425 {
 426   return c->n_obs;
 427 }
 428
 429 double
 430 linreg_sse (const linreg *c)
 431 {
 432   return c->sse;
 433 }
 434
 435 double
 436 linreg_ssreg (const linreg *c)
 437 {
 438   return (c->sst - c->sse);
 439 }
 440
 441 double linreg_sst (const linreg *c)
 442 {
 443   return c->sst;
 444 }
 445
 446 double
 447 linreg_dfmodel ( const linreg *c)
 448 {
 449   return c->dfm;
 450 }
 451
 452 void
 453 linreg_set_depvar_mean (linreg *c, double x)
 454 {
 455   c->depvar_mean = x;
 456 }
 457
 458 double
 459 linreg_get_depvar_mean (linreg *c)
 460 {
 461   return c->depvar_mean;
 462 }