src/math/linreg/linreg.h

   1 /*
   2   lib/linreg/linreg.h
   3
   4   Copyright (C) 2005 Free Software Foundation, Inc. Written by Jason H. Stover.
   5
   6   This program is free software; you can redistribute it and/or modify it under
   7   the terms of the GNU General Public License as published by the Free
   8   Software Foundation; either version 2 of the License, or (at your option)
   9   any later version.
  10
  11   This program is distributed in the hope that it will be useful, but WITHOUT
  12   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13   FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
  14   more details.
  15
  16   You should have received a copy of the GNU General Public License along with
  17   this program; if not, write to the Free Software Foundation, Inc., 51
  18   Franklin Street, Fifth Floor, Boston, MA 02111-1307, USA.
  19 */
  20
  21 #ifndef LINREG_H
  22 #define LINREG_H
  23 #include <stdbool.h>
  24 #include <gsl/gsl_math.h>
  25 #include <gsl/gsl_vector.h>
  26 #include <gsl/gsl_matrix.h>
  27
  28 struct variable;
  29 struct pspp_coeff;
  30 union value;
  31
  32 enum
  33 {
  34   PSPP_LINREG_SWEEP,
  35   PSPP_LINREG_SVD
  36 };
  37
  38
  39
  40 /*
  41   Options describing what special values should be computed.
  42  */
  43 struct pspp_linreg_opts_struct
  44 {
  45   int get_depvar_mean_std;
  46   int *get_indep_mean_std;      /* Array of booleans
  47                                    dictating which
  48                                    independent variables need
  49                                    their means and standard
  50                                    deviations computed within
  51                                    pspp_linreg. This array
  52                                    MUST be of length
  53                                    n_indeps. If element i is
  54                                    1, pspp_linreg will
  55                                    compute the mean and
  56                                    variance of indpendent
  57                                    variable i. If element i
  58                                    is 0, it will not compute
  59                                    the mean and standard
  60                                    deviation, and assume the
  61                                    values are stored.
  62                                    cache->indep_mean[i] is
  63                                    the mean and
  64                                    cache->indep_std[i] is the
  65                                    sample standard deviation. */
  66 };
  67 typedef struct pspp_linreg_opts_struct pspp_linreg_opts;
  68
  69
  70 /*
  71   Find the least-squares estimate of b for the linear model:
  72
  73   Y = Xb + Z
  74
  75   where Y is an n-by-1 column vector, X is an n-by-p matrix of
  76   independent variables, b is a p-by-1 vector of regression coefficients,
  77   and Z is an n-by-1 normally-distributed random vector with independent
  78   identically distributed components with mean 0.
  79
  80   This estimate is found via the sweep operator or singular-value
  81   decomposition with gsl.
  82
  83
  84   References:
  85
  86   1. Matrix Computations, third edition. GH Golub and CF Van Loan.
  87   The Johns Hopkins University Press. 1996. ISBN 0-8018-5414-8.
  88
  89   2. Numerical Analysis for Statisticians. K Lange. Springer. 1999.
  90   ISBN 0-387-94979-8.
  91
  92   3. Numerical Linear Algebra for Applications in Statistics. JE Gentle.
  93   Springer. 1998. ISBN 0-387-98542-5.
  94 */
  95
  96
  97 struct pspp_linreg_cache_struct
  98 {
  99   int n_obs;                    /* Number of observations. */
 100   int n_indeps;                 /* Number of independent variables. */
 101   int n_coeffs;
 102
 103   /*
 104      The variable struct is ignored during estimation. It is here so
 105      the calling procedure can find the variable used in the model.
 106    */
 107   const struct variable *depvar;
 108
 109   gsl_vector *residuals;
 110   struct pspp_coeff **coeff;
 111   int method;                   /* Method to use to estimate parameters. */
 112   /*
 113      Means and standard deviations of the variables.
 114      If these pointers are null when pspp_linreg() is
 115      called, pspp_linreg() will compute their values.
 116
 117      Entry i of indep_means is the mean of independent
 118      variable i, whose observations are stored in the ith
 119      column of the design matrix.
 120    */
 121   double depvar_mean;
 122   double depvar_std;
 123   gsl_vector *indep_means;
 124   gsl_vector *indep_std;
 125
 126   /*
 127      Sums of squares.
 128    */
 129   double ssm;                   /* Sums of squares for the overall model. */
 130   gsl_vector *ss_indeps;        /* Sums of squares from each
 131                                    independent variable. */
 132   double sst;                   /* Sum of squares total. */
 133   double sse;                   /* Sum of squares error. */
 134   double mse;                   /* Mean squared error. This is just sse /
 135                                    dfe, but since it is the best unbiased
 136                                    estimate of the population variance, it
 137                                    has its own entry here. */
 138   gsl_vector *ssx;              /* Centered sums of squares for independent
 139                                    variables, i.e. \sum (x[i] - mean(x))^2. */
 140   double ssy;                   /* Centered sums of squares for dependent
 141                                    variable.
 142                                  */
 143   /*
 144      Covariance matrix of the parameter estimates.
 145    */
 146   gsl_matrix *cov;
 147   /*
 148      Degrees of freedom.
 149    */
 150   double dft;
 151   double dfe;
 152   double dfm;
 153
 154   /*
 155      'Hat' or Hessian matrix, i.e. (X'X)^{-1}, where X is our
 156      design matrix.
 157    */
 158   gsl_matrix *hat;
 159
 160   double (*predict) (const struct variable **, const union value **,
 161                      const void *, int);
 162   double (*residual) (const struct variable **,
 163                       const union value **,
 164                       const union value *, const void *, int);
 165   /*
 166      Returns pointers to the variables used in the model.
 167    */
 168   int (*get_vars) (const void *, const struct variable **);
 169   struct variable *resid;
 170   struct variable *pred;
 171
 172 };
 173
 174 typedef struct pspp_linreg_cache_struct pspp_linreg_cache;
 175
 176
 177
 178 /*
 179   Allocate a pspp_linreg_cache and return a pointer
 180   to it. n is the number of cases, p is the number of
 181   independent variables.
 182  */
 183 pspp_linreg_cache *pspp_linreg_cache_alloc (size_t n, size_t p);
 184
 185 bool pspp_linreg_cache_free (void *);
 186
 187 /*
 188   Fit the linear model via least squares. All pointers passed to pspp_linreg
 189   are assumed to be allocated to the correct size and initialized to the
 190   values as indicated by opts.
 191  */
 192 int
 193 pspp_linreg (const gsl_vector * Y, const gsl_matrix * X,
 194              const pspp_linreg_opts * opts, pspp_linreg_cache * cache);
 195
 196 double
 197 pspp_linreg_predict (const struct variable **, const union value **,
 198                      const void *, int);
 199 double
 200 pspp_linreg_residual (const struct variable **, const union value **,
 201                       const union value *, const void *, int);
 202 /*
 203   All variables used in the model.
 204  */
 205 int pspp_linreg_get_vars (const void *, const struct variable **);
 206 #endif