src/math/linreg.h

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2005 Free Software Foundation, Inc. Written by Jason H. Stover.
   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 #ifndef LINREG_H
  18 #define LINREG_H
  19 #include <stdbool.h>
  20 #include <gsl/gsl_math.h>
  21 #include <gsl/gsl_vector.h>
  22 #include <gsl/gsl_matrix.h>
  23 #include <src/math/coefficient.h>
  24 #include <math/covariance-matrix.h>
  25
  26 enum
  27 {
  28   PSPP_LINREG_CONDITIONAL_INVERSE,
  29   PSPP_LINREG_QR,
  30   PSPP_LINREG_SWEEP,
  31 };
  32
  33
  34
  35 /*
  36   Options describing what special values should be computed.
  37  */
  38 struct pspp_linreg_opts_struct
  39 {
  40   int get_depvar_mean_std;
  41   int *get_indep_mean_std;      /* Array of booleans
  42                                    dictating which
  43                                    independent variables need
  44                                    their means and standard
  45                                    deviations computed within
  46                                    pspp_linreg. This array
  47                                    MUST be of length
  48                                    n_indeps. If element i is
  49                                    1, pspp_linreg will
  50                                    compute the mean and
  51                                    variance of indpendent
  52                                    variable i. If element i
  53                                    is 0, it will not compute
  54                                    the mean and standard
  55                                    deviation, and assume the
  56                                    values are stored.
  57                                    cache->indep_mean[i] is
  58                                    the mean and
  59                                    cache->indep_std[i] is the
  60                                    sample standard deviation. */
  61 };
  62 typedef struct pspp_linreg_opts_struct pspp_linreg_opts;
  63
  64
  65 /*
  66   Find the least-squares estimate of b for the linear model:
  67
  68   Y = Xb + Z
  69
  70   where Y is an n-by-1 column vector, X is an n-by-p matrix of
  71   independent variables, b is a p-by-1 vector of regression coefficients,
  72   and Z is an n-by-1 normally-distributed random vector with independent
  73   identically distributed components with mean 0.
  74
  75   This estimate is found via the sweep operator or singular-value
  76   decomposition with gsl.
  77
  78
  79   References:
  80
  81   1. Matrix Computations, third edition. GH Golub and CF Van Loan.
  82   The Johns Hopkins University Press. 1996. ISBN 0-8018-5414-8.
  83
  84   2. Numerical Analysis for Statisticians. K Lange. Springer. 1999.
  85   ISBN 0-387-94979-8.
  86
  87   3. Numerical Linear Algebra for Applications in Statistics. JE Gentle.
  88   Springer. 1998. ISBN 0-387-98542-5.
  89 */
  90
  91
  92 struct pspp_linreg_cache_struct
  93 {
  94   int n_obs;                    /* Number of observations. */
  95   int n_indeps;                 /* Number of independent variables. */
  96   int n_coeffs;                 /* The intercept is not considered a
  97                                    coefficient here. */
  98
  99   /*
 100     Pointers to the variables.
 101    */
 102   const struct variable *depvar;
 103   const struct variable **indep_vars;
 104
 105   gsl_vector *residuals;
 106   struct pspp_coeff **coeff;
 107   double intercept;
 108   int method;                   /* Method to use to estimate parameters. */
 109   /*
 110      Means and standard deviations of the variables.
 111      If these pointers are null when pspp_linreg() is
 112      called, pspp_linreg() will compute their values.
 113
 114      Entry i of indep_means is the mean of independent
 115      variable i, whose observations are stored in the ith
 116      column of the design matrix.
 117    */
 118   double depvar_mean;
 119   double depvar_std;
 120   gsl_vector *indep_means;
 121   gsl_vector *indep_std;
 122
 123   /*
 124      Sums of squares.
 125    */
 126   double ssm;                   /* Sums of squares for the overall model. */
 127   gsl_vector *ss_indeps;        /* Sums of squares from each
 128                                    independent variable. */
 129   double sst;                   /* Sum of squares total. */
 130   double sse;                   /* Sum of squares error. */
 131   double mse;                   /* Mean squared error. This is just sse /
 132                                    dfe, but since it is the best unbiased
 133                                    estimate of the population variance, it
 134                                    has its own entry here. */
 135   gsl_vector *ssx;              /* Centered sums of squares for independent
 136                                    variables, i.e. \sum (x[i] - mean(x))^2. */
 137   double ssy;                   /* Centered sums of squares for dependent
 138                                    variable.
 139                                  */
 140   /*
 141      Covariance matrix of the parameter estimates.
 142    */
 143   gsl_matrix *cov;
 144   /*
 145      Degrees of freedom.
 146    */
 147   double dft;
 148   double dfe;
 149   double dfm;
 150
 151   /*
 152      'Hat' or Hessian matrix, i.e. (X'X)^{-1}, where X is our
 153      design matrix.
 154    */
 155   gsl_matrix *hat;
 156
 157   double (*predict) (const struct variable **, const union value **,
 158                      const void *, int);
 159   double (*residual) (const struct variable **,
 160                       const union value **,
 161                       const union value *, const void *, int);
 162   /*
 163      Returns pointers to the variables used in the model.
 164    */
 165   int (*get_vars) (const void *, const struct variable **);
 166   struct variable *resid;
 167   struct variable *pred;
 168
 169 };
 170
 171 typedef struct pspp_linreg_cache_struct pspp_linreg_cache;
 172
 173
 174
 175 /*
 176   Allocate a pspp_linreg_cache and return a pointer
 177   to it. n is the number of cases, p is the number of
 178   independent variables.
 179  */
 180 pspp_linreg_cache *pspp_linreg_cache_alloc (const struct variable *, const struct variable **,
 181                                             size_t, size_t);
 182
 183 bool pspp_linreg_cache_free (void *);
 184
 185 /*
 186   Fit the linear model via least squares. All pointers passed to pspp_linreg
 187   are assumed to be allocated to the correct size and initialized to the
 188   values as indicated by opts.
 189  */
 190 int
 191 pspp_linreg (const gsl_vector *, const struct design_matrix *,
 192              const pspp_linreg_opts *, pspp_linreg_cache *);
 193
 194 double
 195 pspp_linreg_predict (const struct variable **, const union value **,
 196                      const void *, int);
 197 double
 198 pspp_linreg_residual (const struct variable **, const union value **,
 199                       const union value *, const void *, int);
 200 /*
 201   All variables used in the model.
 202  */
 203 int pspp_linreg_get_vars (const void *, const struct variable **);
 204
 205 struct pspp_coeff *pspp_linreg_get_coeff (const pspp_linreg_cache
 206                                                        *,
 207                                                        const struct variable
 208                                                        *,
 209                                                        const union value *);
 210 /*
 211   Return or set the standard deviation of the independent variable.
 212  */
 213 double pspp_linreg_get_indep_variable_sd (pspp_linreg_cache *, const struct variable *);
 214 void pspp_linreg_set_indep_variable_sd (pspp_linreg_cache *, const struct variable *, double);
 215 /*
 216   Mean of the independent variable.
 217  */
 218 double pspp_linreg_get_indep_variable_mean (pspp_linreg_cache *, const struct variable *);
 219 void pspp_linreg_set_indep_variable_mean (pspp_linreg_cache *, const struct variable *, double);
 220
 221 /*
 222   Regression using only the covariance matrix.
 223  */
 224 void pspp_linreg_with_cov (const struct covariance_matrix *, pspp_linreg_cache *);
 225 double pspp_linreg_mse (const pspp_linreg_cache *);
 226 #endif