pintos-os.org Git - pspp/blob - lib/linreg/pspp_linreg.h

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