src/math/linreg/linreg.h

   1 /* lib/linreg/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 #ifndef LINREG_H
  23 #define LINREG_H
  24
  25
  26 #include <gsl/gsl_vector.h>
  27 #include <gsl/gsl_matrix.h>
  28
  29 struct variable ;
  30 struct pspp_linreg_coeff;
  31
  32
  33 enum
  34 {
  35   PSPP_LINREG_SWEEP,
  36   PSPP_LINREG_SVD
  37 };
  38
  39
  40
  41 /*
  42   Options describing what special values should be computed.
  43  */
  44 struct pspp_linreg_opts_struct
  45 {
  46   int resid;                    /* Should the residuals be returned? */
  47
  48   int get_depvar_mean_std;
  49   int *get_indep_mean_std;      /* Array of booleans dictating which
  50                                    independent variables need their means
  51                                    and standard deviations computed within
  52                                    pspp_linreg. This array MUST be of
  53                                    length n_indeps. If element i is 1,
  54                                    pspp_linreg will compute the mean and
  55                                    variance of indpendent variable i. If
  56                                    element i is 0, it will not compute the
  57                                    mean and standard deviation, and assume
  58                                    the values are stored.
  59                                    cache->indep_mean[i] is the mean and
  60                                    cache->indep_std[i] is the sample
  61                                    standard deviation.
  62                                  */
  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;
  99
 100   /*
 101      The variable struct is ignored during estimation.
 102      It is here so the calling procedure can
 103      find the variable used in the model.
 104   */
 105   const struct variable *depvar;
 106
 107   gsl_vector *residuals;
 108   struct pspp_linreg_coeff *coeff;
 109   int method;                   /* Method to use to estimate parameters. */
 110   /*
 111      Means and standard deviations of the variables.
 112      If these pointers are null when pspp_linreg() is
 113      called, pspp_linreg() will compute their values.
 114
 115      Entry i of indep_means is the mean of independent
 116      variable i, whose observations are stored in the ith
 117      column of the design matrix.
 118    */
 119   double depvar_mean;
 120   double depvar_std;
 121   gsl_vector *indep_means;
 122   gsl_vector *indep_std;
 123
 124   /*
 125      Sums of squares.
 126    */
 127   double ssm;                   /* Sums of squares for the overall model. */
 128   gsl_vector *ss_indeps;        /* Sums of squares from each
 129                                    independent variable.
 130                                  */
 131   double sst;                   /* Sum of squares total. */
 132   double sse;                   /* Sum of squares error. */
 133   double mse;                   /* Mean squared error. This is just sse / dfe, but
 134                                    since it is the best unbiased estimate of the population
 135                                    variance, it has its own entry here.
 136                                  */
 137   gsl_vector *ssx;              /* Centered sums of squares for independent variables,
 138                                    i.e. \sum (x[i] - mean(x))^2.
 139                                  */
 140   double ssy;                   /* Centered sums of squares for dependent variable. */
 141   /*
 142      Covariance matrix of the parameter estimates.
 143    */
 144   gsl_matrix *cov;
 145   /*
 146      Degrees of freedom.
 147    */
 148   double dft;
 149   double dfe;
 150   double dfm;
 151
 152   /*
 153      'Hat' or Hessian matrix, i.e. (X'X)^{-1}, where X is our
 154      design matrix.
 155    */
 156   gsl_matrix *hat;
 157 };
 158
 159 typedef struct pspp_linreg_cache_struct pspp_linreg_cache;
 160
 161
 162
 163 /*
 164   Allocate a pspp_linreg_cache and return a pointer
 165   to it. n is the number of cases, p is the number of
 166   independent variables.
 167  */
 168 pspp_linreg_cache * pspp_linreg_cache_alloc (size_t n, size_t p);
 169
 170 void pspp_linreg_cache_free (pspp_linreg_cache * c);
 171
 172 /*
 173   Fit the linear model via least squares. All pointers passed to pspp_linreg
 174   are assumed to be allocated to the correct size and initialized to the
 175   values as indicated by opts.
 176  */
 177 int pspp_linreg (const gsl_vector * Y, const gsl_matrix * X,
 178                  const pspp_linreg_opts * opts,
 179                  pspp_linreg_cache * cache);
 180
 181
 182 #endif