a9577d648372dc640056fc07762405ccc0aff6c1
[pspp-builds.git] / 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 #include <src/math/coefficient.h>
24
25 enum
26 {
27   PSPP_LINREG_CONDITIONAL_INVERSE,
28   PSPP_LINREG_QR,
29   PSPP_LINREG_SWEEP,
30 };
31
32
33
34 /*
35   Options describing what special values should be computed.
36  */
37 struct pspp_linreg_opts_struct
38 {
39   int get_depvar_mean_std;
40   int *get_indep_mean_std;      /* Array of booleans
41                                    dictating which
42                                    independent variables need
43                                    their means and standard
44                                    deviations computed within
45                                    pspp_linreg. This array
46                                    MUST be of length
47                                    n_indeps. If element i is
48                                    1, pspp_linreg will
49                                    compute the mean and
50                                    variance of indpendent
51                                    variable i. If element i
52                                    is 0, it will not compute
53                                    the mean and standard
54                                    deviation, and assume the
55                                    values are stored.
56                                    cache->indep_mean[i] is
57                                    the mean and
58                                    cache->indep_std[i] is the
59                                    sample standard deviation. */
60 };
61 typedef struct pspp_linreg_opts_struct pspp_linreg_opts;
62
63
64 /*
65   Find the least-squares estimate of b for the linear model:
66
67   Y = Xb + Z
68
69   where Y is an n-by-1 column vector, X is an n-by-p matrix of
70   independent variables, b is a p-by-1 vector of regression coefficients,
71   and Z is an n-by-1 normally-distributed random vector with independent
72   identically distributed components with mean 0.
73
74   This estimate is found via the sweep operator or singular-value
75   decomposition with gsl.
76
77
78   References:
79
80   1. Matrix Computations, third edition. GH Golub and CF Van Loan.
81   The Johns Hopkins University Press. 1996. ISBN 0-8018-5414-8.
82
83   2. Numerical Analysis for Statisticians. K Lange. Springer. 1999.
84   ISBN 0-387-94979-8.
85
86   3. Numerical Linear Algebra for Applications in Statistics. JE Gentle.
87   Springer. 1998. ISBN 0-387-98542-5.
88 */
89
90
91 struct pspp_linreg_cache_struct
92 {
93   int n_obs;                    /* Number of observations. */
94   int n_indeps;                 /* Number of independent variables. */
95   int n_coeffs;                 /* The intercept is not considered a
96                                    coefficient here. */
97
98   /*
99     Pointers to the variables.
100    */
101   const struct variable *depvar;
102   const struct variable **indep_vars;
103
104   gsl_vector *residuals;
105   struct pspp_coeff **coeff;
106   double intercept;
107   int method;                   /* Method to use to estimate parameters. */
108   /*
109      Means and standard deviations of the variables.
110      If these pointers are null when pspp_linreg() is
111      called, pspp_linreg() will compute their values.
112
113      Entry i of indep_means is the mean of independent
114      variable i, whose observations are stored in the ith
115      column of the design matrix.
116    */
117   double depvar_mean;
118   double depvar_std;
119   gsl_vector *indep_means;
120   gsl_vector *indep_std;
121
122   /*
123      Sums of squares.
124    */
125   double ssm;                   /* Sums of squares for the overall model. */
126   gsl_vector *ss_indeps;        /* Sums of squares from each
127                                    independent variable. */
128   double sst;                   /* Sum of squares total. */
129   double sse;                   /* Sum of squares error. */
130   double mse;                   /* Mean squared error. This is just sse /
131                                    dfe, but since it is the best unbiased
132                                    estimate of the population variance, it
133                                    has its own entry here. */
134   gsl_vector *ssx;              /* Centered sums of squares for independent
135                                    variables, i.e. \sum (x[i] - mean(x))^2. */
136   double ssy;                   /* Centered sums of squares for dependent
137                                    variable.
138                                  */
139   /*
140      Covariance matrix of the parameter estimates.
141    */
142   gsl_matrix *cov;
143   /*
144      Degrees of freedom.
145    */
146   double dft;
147   double dfe;
148   double dfm;
149
150   /*
151      'Hat' or Hessian matrix, i.e. (X'X)^{-1}, where X is our
152      design matrix.
153    */
154   gsl_matrix *hat;
155
156   double (*predict) (const struct variable **, const union value **,
157                      const void *, int);
158   double (*residual) (const struct variable **,
159                       const union value **,
160                       const union value *, const void *, int);
161   /*
162      Returns pointers to the variables used in the model.
163    */
164   int (*get_vars) (const void *, const struct variable **);
165   struct variable *resid;
166   struct variable *pred;
167
168 };
169
170 typedef struct pspp_linreg_cache_struct pspp_linreg_cache;
171
172
173
174 /*
175   Allocate a pspp_linreg_cache and return a pointer
176   to it. n is the number of cases, p is the number of
177   independent variables.
178  */
179 pspp_linreg_cache *pspp_linreg_cache_alloc (const struct variable *, const struct variable **, 
180                                             size_t, size_t);
181
182 bool pspp_linreg_cache_free (void *);
183
184 /*
185   Fit the linear model via least squares. All pointers passed to pspp_linreg
186   are assumed to be allocated to the correct size and initialized to the
187   values as indicated by opts.
188  */
189 int
190 pspp_linreg (const gsl_vector *, const struct design_matrix *,
191              const pspp_linreg_opts *, pspp_linreg_cache *);
192
193 double
194 pspp_linreg_predict (const struct variable **, const union value **,
195                      const void *, int);
196 double
197 pspp_linreg_residual (const struct variable **, const union value **,
198                       const union value *, const void *, int);
199 /*
200   All variables used in the model.
201  */
202 int pspp_linreg_get_vars (const void *, struct variable **);
203
204 struct pspp_coeff *pspp_linreg_get_coeff (const pspp_linreg_cache
205                                                        *,
206                                                        const struct variable
207                                                        *,
208                                                        const union value *);
209 /*
210   Return or set the standard deviation of the independent variable.
211  */
212 double pspp_linreg_get_indep_variable_sd (pspp_linreg_cache *, const struct variable *);
213 void pspp_linreg_set_indep_variable_sd (pspp_linreg_cache *, const struct variable *, double);
214 /*
215   Mean of the independent variable.
216  */
217 double pspp_linreg_get_indep_variable_mean (pspp_linreg_cache *, const struct variable *);
218 void pspp_linreg_set_indep_variable_mean (pspp_linreg_cache *, const struct variable *, double);
219
220 /*
221   Regression using only the covariance matrix.
222  */
223 int pspp_linreg_with_cov (const struct design_matrix *, pspp_linreg_cache *);
224 #endif