src/math/linreg.h

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2005 Free Software Foundation, Inc. Written by Jason H. Stover.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 #ifndef LINREG_H
  18 #define LINREG_H
  19 #include <stdbool.h>
  20 #include <gsl/gsl_math.h>
  21 #include <gsl/gsl_vector.h>
  22 #include <gsl/gsl_matrix.h>
  23
  24 enum
  25 {
  26   LINREG_CONDITIONAL_INVERSE,
  27   LINREG_QR,
  28   LINREG_SWEEP,
  29 };
  30
  31
  32
  33 /*
  34   Options describing what special values should be computed.
  35  */
  36 struct pspp_linreg_opts_struct
  37 {
  38   int get_depvar_mean_std;
  39   int *get_indep_mean_std;      /* Array of booleans
  40                                    dictating which
  41                                    independent variables need
  42                                    their means and standard
  43                                    deviations computed within
  44                                    pspp_linreg. This array
  45                                    MUST be of length
  46                                    n_indeps. If element i is
  47                                    1, pspp_linreg will
  48                                    compute the mean and
  49                                    variance of indpendent
  50                                    variable i. If element i
  51                                    is 0, it will not compute
  52                                    the mean and standard
  53                                    deviation, and assume the
  54                                    values are stored.
  55                                    cache->indep_mean[i] is
  56                                    the mean and
  57                                    cache->indep_std[i] is the
  58                                    sample standard deviation. */
  59 };
  60 typedef struct pspp_linreg_opts_struct pspp_linreg_opts;
  61
  62
  63 /*
  64   Find the least-squares estimate of b for the linear model:
  65
  66   Y = Xb + Z
  67
  68   where Y is an n-by-1 column vector, X is an n-by-p matrix of
  69   independent variables, b is a p-by-1 vector of regression coefficients,
  70   and Z is an n-by-1 normally-distributed random vector with independent
  71   identically distributed components with mean 0.
  72
  73   This estimate is found via the sweep operator or singular-value
  74   decomposition with gsl.
  75
  76
  77   References:
  78
  79   1. Matrix Computations, third edition. GH Golub and CF Van Loan.
  80   The Johns Hopkins University Press. 1996. ISBN 0-8018-5414-8.
  81
  82   2. Numerical Analysis for Statisticians. K Lange. Springer. 1999.
  83   ISBN 0-387-94979-8.
  84
  85   3. Numerical Linear Algebra for Applications in Statistics. JE Gentle.
  86   Springer. 1998. ISBN 0-387-98542-5.
  87 */
  88
  89
  90 struct linreg_struct
  91 {
  92   double n_obs;                 /* Number of observations. */
  93   int n_indeps;                 /* Number of independent variables. */
  94   int n_coeffs;                 /* The intercept is not considered a
  95                                    coefficient here. */
  96
  97   /*
  98     Pointers to the variables.
  99    */
 100   const struct variable *depvar;
 101   const struct variable **indep_vars;
 102
 103   double *coeff;
 104   double intercept;
 105   int method;                   /* Method to use to estimate parameters. */
 106   /*
 107      Means and standard deviations of the variables.
 108      If these pointers are null when pspp_linreg() is
 109      called, pspp_linreg() will compute their values.
 110
 111      Entry i of indep_means is the mean of independent
 112      variable i, whose observations are stored in the ith
 113      column of the design matrix.
 114    */
 115   double depvar_mean;
 116   double depvar_std;
 117   gsl_vector *indep_means;
 118   gsl_vector *indep_std;
 119
 120   /*
 121      Sums of squares.
 122    */
 123   double ssm;                   /* Sums of squares for the overall model. */
 124   gsl_vector *ss_indeps;        /* Sums of squares from each
 125                                    independent variable. */
 126   double sst;                   /* Sum of squares total. */
 127   double sse;                   /* Sum of squares error. */
 128   double mse;                   /* Mean squared error. This is just sse /
 129                                    dfe, but since it is the best unbiased
 130                                    estimate of the population variance, it
 131                                    has its own entry here. */
 132   gsl_vector *ssx;              /* Centered sums of squares for independent
 133                                    variables, i.e. \sum (x[i] - mean(x))^2. */
 134   double ssy;                   /* Centered sums of squares for dependent
 135                                    variable.
 136                                  */
 137   /*
 138      Covariance matrix of the parameter estimates.
 139    */
 140   gsl_matrix *cov;
 141   /*
 142      Degrees of freedom.
 143    */
 144   double dft;
 145   double dfe;
 146   double dfm;
 147
 148   /*
 149      'Hat' or Hessian matrix, i.e. (X'X)^{-1}, where X is our
 150      design matrix.
 151    */
 152   gsl_matrix *hat;
 153
 154   struct variable *pred;
 155   struct variable *resid;
 156 };
 157
 158 typedef struct linreg_struct linreg;
 159
 160
 161
 162 linreg *linreg_alloc (const struct variable *, const struct variable **,
 163                       double, size_t);
 164
 165 bool linreg_free (void *);
 166
 167 /*
 168   Fit the linear model via least squares. All pointers passed to pspp_linreg
 169   are assumed to be allocated to the correct size and initialized to the
 170   values as indicated by opts.
 171  */
 172 void
 173 linreg_fit (const gsl_matrix *, linreg *);
 174
 175 double
 176 linreg_predict (const linreg *, const double *, size_t);
 177 double
 178 linreg_residual (const linreg *, double, const double *, size_t);
 179 const struct variable ** linreg_get_vars (const linreg *);
 180
 181 /*
 182   Return or set the standard deviation of the independent variable.
 183  */
 184 double linreg_get_indep_variable_sd (linreg *, size_t);
 185 void linreg_set_indep_variable_sd (linreg *, size_t, double);
 186 /*
 187   Mean of the independent variable.
 188  */
 189 double linreg_get_indep_variable_mean (linreg *, size_t);
 190 void linreg_set_indep_variable_mean (linreg *, size_t, double);
 191
 192 double linreg_mse (const linreg *);
 193
 194 double linreg_intercept (const linreg *);
 195
 196 gsl_matrix * linreg_cov (const linreg *);
 197 double linreg_coeff (const linreg *, size_t);
 198 const struct variable * linreg_indep_var (const linreg *, size_t);
 199 size_t linreg_n_coeffs (const linreg *);
 200 size_t linreg_n_obs (const linreg *);
 201 double linreg_sse (const linreg *);
 202 double linreg_ssreg (const linreg *);
 203 double linreg_dfmodel (const linreg *);
 204 double linreg_sst (const linreg *);
 205 #endif