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