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 enum
  58 {
  59   PSPP_LINREG_SWEEP,
  60   PSPP_LINREG_SVD
  61 };
  62
  63 /*
  64   Cache for the relevant data from the model. There are several
  65   members which the caller may not use, and which could use a lot of
  66   storage. Therefore non-essential members of the struct will be
  67   allocated only when requested.
  68  */
  69 struct pspp_linreg_cache_struct
  70 {
  71   int n_obs;                    /* Number of observations. */
  72   int n_indeps;                 /* Number of independent variables. */
  73   gsl_vector *depvar;
  74   gsl_matrix *indepvar;
  75   gsl_vector *residuals;
  76   gsl_vector *param_estimates;
  77   int method;                   /* Method to use to estimate parameters. */
  78   /*
  79      Means and standard deviations of the variables.
  80      If these pointers are null when pspp_linreg() is
  81      called, pspp_linreg() will compute their values.
  82
  83      Entry i of indep_means is the mean of independent
  84      variable i, whose observations are stored in column i
  85      of indepvar.
  86    */
  87   double depvar_mean;
  88   double depvar_std;
  89   gsl_vector *indep_means;
  90   gsl_vector *indep_std;
  91
  92   /*
  93      Sums of squares.
  94    */
  95   double ssm;                   /* Sums of squares for the overall model. */
  96   gsl_vector *ss_indeps;        /* Sums of squares from each
  97                                    independent variable.
  98                                  */
  99   double sst;                   /* Sum of squares total. */
 100   double sse;                   /* Sum of squares error. */
 101   double mse;                   /* Mean squared error. This is just sse / dfe, but
 102                                    since it is the best unbiased estimate of the population
 103                                    variance, it has its own entry here.
 104                                  */
 105   gsl_vector *ssx;              /* Centered sums of squares for independent variables,
 106                                    i.e. \sum (x[i] - mean(x))^2.
 107                                  */
 108   double ssy;                   /* Centered sums of squares for dependent variable. */
 109   /*
 110      Covariance matrix of the parameter estimates.
 111    */
 112   gsl_matrix *cov;
 113   /*
 114      Degrees of freedom.
 115    */
 116   double dft;
 117   double dfe;
 118   double dfm;
 119
 120   /*
 121      'Hat' or Hessian matrix, i.e. (X'X)^{-1}, where X is our
 122      design matrix.
 123    */
 124   gsl_matrix *hat;
 125 };
 126 typedef struct pspp_linreg_cache_struct pspp_linreg_cache;
 127
 128 /*
 129   Options describing what special values should be computed.
 130  */
 131 struct pspp_linreg_opts_struct
 132 {
 133   int resid;                    /* Should the residuals be returned? */
 134
 135   int get_depvar_mean_std;
 136   int *get_indep_mean_std;      /* Array of booleans dictating which
 137                                    independent variables need their means
 138                                    and standard deviations computed within
 139                                    pspp_linreg. This array MUST be of
 140                                    length n_indeps. If element i is 1,
 141                                    pspp_linreg will compute the mean and
 142                                    variance of indpendent variable i. If
 143                                    element i is 0, it will not compute the
 144                                    mean and standard deviation, and assume
 145                                    the values are stored.
 146                                    cache->indep_mean[i] is the mean and
 147                                    cache->indep_std[i] is the sample
 148                                    standard deviation.
 149                                  */
 150 };
 151 typedef struct pspp_linreg_opts_struct pspp_linreg_opts;
 152
 153 int pspp_reg_sweep (gsl_matrix * A);
 154
 155 pspp_linreg_cache *pspp_linreg_cache_alloc (size_t n, size_t p);
 156
 157 void pspp_linreg_cache_free (pspp_linreg_cache * cache);
 158
 159 int pspp_linreg (const gsl_vector * Y, const gsl_matrix * X,
 160                  const pspp_linreg_opts * opts, pspp_linreg_cache * cache);
 161 #endif