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