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

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2005, 2010, 2011, 2017 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, bool origin)
  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;
  95   if (!origin)
  96     c->dft--;
  97
  98   c->dfm = p;
  99   c->dfe = c->dft - c->dfm;
 100   c->intercept = 0.0;
 101   c->depvar_mean = 0.0;
 102   /*
 103      Default settings.
 104    */
 105   c->method = LINREG_SWEEP;
 106
 107   c->refcnt = 1;
 108
 109   c->origin = origin;
 110
 111   return c;
 112 }
 113
 114
 115 void
 116 linreg_ref (linreg *c)
 117 {
 118   c->refcnt++;
 119 }
 120
 121 void
 122 linreg_unref (linreg *c)
 123 {
 124   if (--c->refcnt == 0)
 125     {
 126       gsl_vector_free (c->indep_means);
 127       gsl_vector_free (c->indep_std);
 128       gsl_matrix_free (c->cov);
 129       free (c->indep_vars);
 130       free (c->coeff);
 131       free (c);
 132     }
 133 }
 134
 135 static void
 136 post_sweep_computations (linreg *l, gsl_matrix *sw)
 137 {
 138   double m;
 139   size_t i;
 140   size_t j;
 141   int rc;
 142
 143   assert (sw != NULL);
 144   assert (l != NULL);
 145
 146   l->sse = gsl_matrix_get (sw, l->n_indeps, l->n_indeps);
 147   l->mse = l->sse / l->dfe;
 148   /*
 149     Get the intercept.
 150   */
 151   m = l->depvar_mean;
 152   for (i = 0; i < l->n_indeps; i++)
 153     {
 154       double tmp = gsl_matrix_get (sw, i, l->n_indeps);
 155       l->coeff[i] = tmp;
 156       m -= tmp * linreg_get_indep_variable_mean (l, i);
 157     }
 158   /*
 159     Get the covariance matrix of the parameter estimates.
 160     Only the upper triangle is necessary.
 161   */
 162
 163   /*
 164     The loops below do not compute the entries related
 165     to the estimated intercept.
 166   */
 167   for (i = 0; i < l->n_indeps; i++)
 168     for (j = i; j < l->n_indeps; j++)
 169       {
 170         double tmp = -1.0 * l->mse * gsl_matrix_get (sw, i, j);
 171         gsl_matrix_set (l->cov, i + 1, j + 1, tmp);
 172       }
 173
 174   if (! l->origin)
 175     {
 176       gsl_matrix *xm;
 177       gsl_matrix_view xtx;
 178       gsl_matrix_view xmxtx;
 179       /*
 180         Get the covariances related to the intercept.
 181       */
 182       xtx = gsl_matrix_submatrix (sw, 0, 0, l->n_indeps, l->n_indeps);
 183       xmxtx = gsl_matrix_submatrix (l->cov, 0, 1, 1, l->n_indeps);
 184       xm = gsl_matrix_calloc (1, l->n_indeps);
 185       for (i = 0; i < xm->size2; i++)
 186         {
 187           gsl_matrix_set (xm, 0, i,
 188                           linreg_get_indep_variable_mean (l, i));
 189         }
 190       rc = gsl_blas_dsymm (CblasRight, CblasUpper, l->mse,
 191                            &xtx.matrix, xm, 0.0, &xmxtx.matrix);
 192       gsl_matrix_free (xm);
 193       if (rc == GSL_SUCCESS)
 194         {
 195           double tmp = l->mse / l->n_obs;
 196           for (i = 1; i < 1 + l->n_indeps; i++)
 197             {
 198               tmp -= gsl_matrix_get (l->cov, 0, i)
 199                 * linreg_get_indep_variable_mean (l, i - 1);
 200             }
 201           gsl_matrix_set (l->cov, 0, 0, tmp);
 202
 203           l->intercept = m;
 204         }
 205       else
 206         {
 207           fprintf (stderr, "%s:%d:gsl_blas_dsymm: %s\n",
 208                    __FILE__, __LINE__, gsl_strerror (rc));
 209           exit (rc);
 210         }
 211     }
 212 }
 213
 214 /*
 215   Predict the value of the dependent variable with the new set of
 216   predictors. VALS are assumed to be in the order corresponding to the
 217   order of the coefficients in the linreg struct.
 218  */
 219 double
 220 linreg_predict (const linreg *c, const double *vals, size_t n_vals)
 221 {
 222   size_t j;
 223   double result;
 224
 225   assert (n_vals = c->n_coeffs);
 226   if (vals == NULL || c == NULL)
 227     {
 228       return GSL_NAN;
 229     }
 230   if (c->coeff == NULL)
 231     {
 232       /* The stupid model: just guess the mean. */
 233       return c->depvar_mean;
 234     }
 235   result = c->intercept;
 236
 237   for (j = 0; j < n_vals; j++)
 238     {
 239       result += linreg_coeff (c, j) * vals[j];
 240     }
 241
 242   return result;
 243 }
 244
 245 double
 246 linreg_residual (const linreg *c, double obs, const double *vals, size_t n_vals)
 247 {
 248   if (vals == NULL || c == NULL)
 249     {
 250       return GSL_NAN;
 251     }
 252   return (obs - linreg_predict (c, vals, n_vals));
 253 }
 254
 255 /*
 256   Mean of the independent variable.
 257  */
 258 double linreg_get_indep_variable_mean (const 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   double intcpt_coef = 0.0;
 274   double intercept_variance = 0.0;
 275   gsl_matrix *xtx;
 276   gsl_matrix *q;
 277   gsl_matrix *r;
 278   gsl_vector *xty;
 279   gsl_vector *tau;
 280   gsl_vector *params;
 281   size_t i;
 282   size_t j;
 283
 284   xtx = gsl_matrix_alloc (cov->size1 - 1, cov->size2 - 1);
 285   xty = gsl_vector_alloc (cov->size1 - 1);
 286   tau = gsl_vector_alloc (cov->size1 - 1);
 287   params = gsl_vector_alloc (cov->size1 - 1);
 288
 289   for (i = 0; i < xtx->size1; i++)
 290     {
 291       gsl_vector_set (xty, i, gsl_matrix_get (cov, cov->size2 - 1, i));
 292       for (j = 0; j < xtx->size2; j++)
 293         {
 294           gsl_matrix_set (xtx, i, j, gsl_matrix_get (cov, i, j));
 295         }
 296     }
 297   gsl_linalg_QR_decomp (xtx, tau);
 298   q = gsl_matrix_alloc (xtx->size1, xtx->size2);
 299   r = gsl_matrix_alloc (xtx->size1, xtx->size2);
 300
 301   gsl_linalg_QR_unpack (xtx, tau, q, r);
 302   gsl_linalg_QR_solve (xtx, tau, xty, params);
 303   for (i = 0; i < params->size; i++)
 304     {
 305       l->coeff[i] = gsl_vector_get (params, i);
 306     }
 307   l->sst = gsl_matrix_get (cov, cov->size1 - 1, cov->size2 - 1);
 308   l->ssm = 0.0;
 309   for (i = 0; i < l->n_indeps; i++)
 310     {
 311       l->ssm += gsl_vector_get (xty, i) * l->coeff[i];
 312     }
 313   l->sse = l->sst - l->ssm;
 314
 315   gsl_blas_dtrsm (CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit, linreg_mse (l),
 316                   r, q);
 317   /* Copy the lower triangle into the upper triangle. */
 318   for (i = 0; i < q->size1; i++)
 319     {
 320       gsl_matrix_set (l->cov, i + 1, i + 1, gsl_matrix_get (q, i, i));
 321       for (j = i + 1; j < q->size2; j++)
 322         {
 323           intercept_variance -= 2.0 * gsl_matrix_get (q, i, j) *
 324             linreg_get_indep_variable_mean (l, i) *
 325             linreg_get_indep_variable_mean (l, j);
 326           gsl_matrix_set (q, i, j, gsl_matrix_get (q, j, i));
 327         }
 328     }
 329   l->intercept = linreg_get_depvar_mean (l);
 330   for (i = 0; i < l->n_indeps; i++)
 331     {
 332       double 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, l->dependent_column);
 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 const 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 (const linreg *c)
 460 {
 461   return c->depvar_mean;
 462 }