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