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