linear-least-squares.lisp 15.8 KB
Newer Older
1
2
;; Linear least squares, or linear regression
;; Liam Healy <2008-01-21 12:41:46EST linear-least-squares.lisp>
Liam M. Healy's avatar
Liam M. Healy committed
3
;; Time-stamp: <2012-12-31 12:30:36EST linear-least-squares.lisp>
4
;;
Liam M. Healy's avatar
Liam M. Healy committed
5
;; Copyright 2008, 2009, 2010, 2011, 2012 Liam M. Healy
6
7
8
9
10
11
12
13
14
15
16
17
18
19
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <http://www.gnu.org/licenses/>.
20
21
22

(in-package :gsl)

23
;;; /usr/include/gsl/gsl_fit.h
24
;;; /usr/include/gsl/gsl_multifit.h
25

26
27
28
29
30
;;;;****************************************************************************
;;;; Linear regression
;;;;****************************************************************************

;;; Error in GSL documentation or in code: vectors must be equal
Liam M. Healy's avatar
Liam M. Healy committed
31
;;; lengths but strides are not?  Shouldn't (floor length stride) be n?
32
33
34
35
;;; Error in GSL documentation for gsl_fit_linear_est, "c00" instead
;;; of "cov00" etc.  Last arg to gsl_fit_wmul is labelled sumsq but
;;; referred to as chisq.

36
37
38
(defmfun linear-fit
    (x y &optional weight (x-stride 1) (y-stride 1) (weight-stride 1))
  ("gsl_fit_linear" "gsl_fit_wlinear")
39
40
41
  ((((grid:foreign-pointer x) :pointer) (x-stride :sizet)
    ((grid:foreign-pointer y) :pointer) (y-stride :sizet)
    ((dim0 x) :sizet)
42
43
44
    (c0 (:pointer :double)) (c1 (:pointer :double))
    (cov00 (:pointer :double)) (cov01 (:pointer :double))
    (cov11 (:pointer :double)) (sumsq (:pointer :double)))
45
46
47
48
   (((grid:foreign-pointer x) :pointer) (x-stride :sizet)
    ((grid:foreign-pointer weight) :pointer) (weight-stride :sizet)
    ((grid:foreign-pointer y) :pointer) (y-stride :sizet)
    ((dim0 x) :sizet)
49
50
51
    (c0 (:pointer :double)) (c1 (:pointer :double))
    (cov00 (:pointer :double)) (cov01 (:pointer :double))
    (cov11 (:pointer :double)) (sumsq (:pointer :double))))
52
  :inputs (x weight y)
53
  :switch (weight weight-stride)
54
55
  :documentation			; FDL
  "Compute the best-fit linear regression coefficients
56
   c0, c1 of the model Y = c_0 + c_1 X for the weighted or unweighted
57
58
   dataset (x, y), two vectors of equal length with strides
   x-stride and y-stride, and return as the first two values.
59
   The vector weight if given, of the same length
60
61
62
63
64
   and stride w-stride, specifies the weight of each datapoint. The
   weight is the reciprocal of the variance for each datapoint in y.

   The covariance matrix for the parameters (c0, c1) is
   computed using the weights and returned via the parameters
65
66
67
   (cov00, cov01, c0v01) as the next three values.  The weighted or
   unweighted sum of squares of the residuals from the best-fit line,
   \chi^2, is returned in as the last value.")
68

liam's avatar
liam committed
69
(defmfun linear-estimate (x c0 c1 cov00 cov01 cov11)
70
71
72
  "gsl_fit_linear_est"
  ((x :double) (c0 :double) (c1 :double)
   (cov00 :double) (cov01 :double) (cov11 :double)
73
   (y (:pointer :double)) (y-error (:pointer :double)))
74
75
76
77
78
79
80
81
82
83
84
  :documentation			; FDL
  "Use the best-fit linear regression coefficients
   c0, c1 and their covariance
   cov00, cov01, cov11 to compute the fitted function
   y and its standard deviation y-error for the model
   Y = c_0 + c_1 X at the point x.")

;;;;****************************************************************************
;;;; Linear fitting without a constant term
;;;;****************************************************************************

85
86
87
(defmfun multiplier-fit
    (x y &optional weight (x-stride 1) (y-stride 1) (weight-stride 1))
  ("gsl_fit_mul" "gsl_fit_wmul")
88
89
90
  ((((grid:foreign-pointer x) :pointer) (x-stride :sizet)
    ((grid:foreign-pointer y) :pointer) (y-stride :sizet)
    ((dim0 x) :sizet)
91
92
    (c1 (:pointer :double)) (cov11 (:pointer :double))
    (sumsq (:pointer :double)))
93
94
95
96
   (((grid:foreign-pointer x) :pointer) (x-stride :sizet)
    ((grid:foreign-pointer weight) :pointer) (weight-stride :sizet)
    ((grid:foreign-pointer y) :pointer) (y-stride :sizet)
    ((dim0 x) :sizet)
97
98
    (c1 (:pointer :double)) (cov11 (:pointer :double))
    (sumsq (:pointer :double))))
99
100
  :inputs (x y weight)
  :switch (weight weight-stride)
101
102
  :documentation			; FDL
  "Compute the best-fit linear regression coefficient
103
   c1 of the model Y = c_1 X for the weighted or unweighted datasets
104
105
106
107
108
109
110
111
112
113
   (x, y), two vectors of equal length with strides
   x-stride and y-stride.  The vector weight of the same length
   and of stride w-stide specifies the weight of each datapoint. The
   weight is the reciprocal of the variance for each datapoint in y.
 
   The variance of the parameter c1 is computed using the weights
   and returned as the second value.  The weighted sum of
   squares of the residuals from the best-fit line, \chi^2, is
   returned as the last value.")

liam's avatar
liam committed
114
(defmfun multiplier-estimate (x c1 cov11)
115
116
  "gsl_fit_mul_est"
  ((x :double) (c1 :double) (cov11 :double)
117
   (y (:pointer :double)) (y-error (:pointer :double)))
118
119
120
121
122
123
124
125
126
127
  :documentation			; FDL
  "Use the best-fit linear regression coefficient
   c1 and its covariance cov11 to compute the fitted function
   y and its standard deviation y-error for the model
   Y = c_0 + c_1 X at the point x.")

;;;;****************************************************************************
;;;; Multiparameter fitting
;;;;****************************************************************************

128
129
(defmobject fit-workspace
    "gsl_multifit_linear"
130
  ((number-of-observations :sizet) (number-of-parameters :sizet))
131
  "multi-dimensional root solver with function only"
Liam Healy's avatar
Liam Healy committed
132
  :documentation
133
  "Make a workspace for a multidimensional linear least-squares fit.")
134

135
136
137
138
139
140
(defun size-array (array-or-size)
  (if (numberp array-or-size)
      array-or-size
      (dim0 array-or-size)))

(defun default-covariance (parameters-or-size)
141
  (grid:make-foreign-array
142
143
144
145
146
147
148
149
   'double-float
   :dimensions
   (let ((s (size-array parameters-or-size))) (list s s))))

(defun default-lls-workspace (observations parameters-or-size)
  (make-fit-workspace
   (dim0 observations) (size-array parameters-or-size)))

150
151
152
(export 'linear-mfit)
(defun linear-mfit
    (model observations parameters-or-size
153
     &optional weight tolerance
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
     (covariance (default-covariance parameters-or-size))
     (workspace (default-lls-workspace observations parameters-or-size)))
  "Compute the best-fit parameters c of the weighted or unweighted
   model y = X c for the observations y and optional weights
   and the model matrix X.  The covariance matrix of
   the model parameters is computed with the given weights.  The
   weighted sum of squares of the residuals from the best-fit,
   chi^2, is returned as the last value.

   The best-fit is found by singular value decomposition of the matrix
   model using the preallocated workspace provided. The modified
   Golub-Reinsch SVD algorithm is used for the unweighted solution,
   with column scaling to improve the accuracy of the singular values.
   Any components which have zero singular value (to machine
   precision) are discarded from the fit.

   If tolerance is a double-float, the SVD algorithm is used.
   If it is nil the non-svd algorithm is used."
  (if tolerance
      (linear-mfit-svd
       model observations parameters-or-size
       tolerance weight covariance workspace)
      (linear-mfit-nosvd
       model observations parameters-or-size
       weight covariance workspace)))

(defmfun linear-mfit-nosvd
181
    (model observations parameters-or-size
182
	   &optional
183
	   weight
184
185
	   (covariance (default-covariance parameters-or-size))
	   (workspace (default-lls-workspace observations parameters-or-size))
186
187
	   &aux
	   (parameters (vdf parameters-or-size)))
188
  ("gsl_multifit_linear" "gsl_multifit_wlinear")
189
190
191
192
  ((((mpointer model) :pointer)
    ((mpointer observations) :pointer)
    ((mpointer parameters) :pointer)
    ((mpointer covariance) :pointer)
193
    (chisq (:pointer :double))
194
195
196
197
198
    ((mpointer workspace) :pointer))
   (((mpointer model) :pointer)
    ((mpointer weight) :pointer)
    ((mpointer observations) :pointer)
    ((mpointer parameters) :pointer)
199
    ((mpointer covariance) :pointer)
200
    (chisq (:pointer :double))
201
    ((mpointer workspace) :pointer)))
202
  :inputs (model weight observations)
203
  :outputs (parameters covariance)
204
  :switch (weight)
205
  :return (parameters covariance (cffi:mem-ref chisq :double))
206
207
  :export nil
  :index linear-mfit
208
  :documentation			; FDL
209
210
  "Compute the best-fit parameters c of the weighted or unweighted
   model y = X c for the observations y and optional weights
211
212
213
   and the model matrix X.  The covariance matrix of
   the model parameters is computed with the given weights.  The
   weighted sum of squares of the residuals from the best-fit,
liam's avatar
liam committed
214
   chi^2, is returned as the last value.
215
216

   The best-fit is found by singular value decomposition of the matrix
217
218
219
220
221
   model using the preallocated workspace provided. The modified
   Golub-Reinsch SVD algorithm is used for the unweighted solution,
   with column scaling to improve the accuracy of the singular values.
   Any components which have zero singular value (to machine
   precision) are discarded from the fit.")
222

223
224
225
(defmfun linear-mfit-svd
    (model observations parameters-or-size tolerance
	   &optional weight
226
227
228
	   (covariance (default-covariance parameters-or-size))
	   (workspace (default-lls-workspace observations parameters-or-size))
	   &aux (parameters (vdf parameters-or-size)))
229
  ("gsl_multifit_linear_svd" "gsl_multifit_wlinear_svd")
230
231
  ((((mpointer model) :pointer)
    ((mpointer observations) :pointer)
232
    (tolerance :double)
233
    (rank (:pointer :sizet))
234
235
    ((mpointer parameters) :pointer)
    ((mpointer covariance) :pointer)
236
    (chisq (:pointer :double))
237
238
239
240
241
    ((mpointer workspace) :pointer))
   (((mpointer model) :pointer)
    ((mpointer weight) :pointer)
    ((mpointer observations) :pointer)
    (tolerance :double)
242
    (rank (:pointer :sizet))
243
244
    ((mpointer parameters) :pointer)
    ((mpointer covariance) :pointer)
245
    (chisq (:pointer :double))
246
    ((mpointer workspace) :pointer)))
247
248
  :inputs (model weight observations)
  :outputs (parameters covariance)
249
  :switch (weight)
250
  :return ((cffi:mem-ref chisq :double) (cffi:mem-ref rank ':sizet))
251
252
  :export nil
  :index linear-mfit
253
  :documentation			; FDL
254
255
256
257
258
259
  "Compute the best-fit parameters c of the weighted or unweighted
   model y = X c for the observations y and weights and the model
   matrix X.  The covariance matrix of the model parameters is
   computed with the given weights.  The weighted or unweighted sum of
   squares of the residuals from the best-fit, chi^2, is returned as
   the first value.
260
261

   The best-fit is found by singular value decomposition of the matrix
262
263
264
265
266
267
268
269
   model using the preallocated workspace provided.  The modified
   Golub-Reinsch SVD algorithm is used, with column scaling to improve
   the accuracy of the singular values (unweighted).  Any components
   which have zero singular value (to machine precision) are discarded
   from the fit.  In the second form of the function the components
   are discarded if the ratio of singular values s_i/s_0 falls below
   the user-specified tolerance, and the effective rank is returned as
   the second value.")
270

liam's avatar
liam committed
271
(defmfun multi-linear-estimate (x coefficients covariance)
272
  "gsl_multifit_linear_est"
273
  (((mpointer x) :pointer) ((mpointer coefficients) :pointer)
274
275
   ((mpointer covariance) :pointer)
   (y (:pointer :double)) (y-error (:pointer :double)))
276
  :inputs (x coefficients covariance)
277
278
279
280
281
282
  :documentation			; FDL
  "Use the best-fit multilinear regression coefficients
   and their covariance matrix to compute the fitted function value
   y and its standard deviation for the model y = x.c
   at the point x.")

283
284
285
286
(defmfun multi-linear-residuals
    (x observations coefficients
       &optional
       (residuals
287
	(grid:make-foreign-array 'double-float :dimensions (grid:dimensions observations))))
288
289
290
291
292
293
294
295
296
297
  "gsl_multifit_linear_residuals"
  (((mpointer x) :pointer) ((mpointer observations) :pointer)
   ((mpointer coefficients) :pointer) ((mpointer residuals) :pointer))
  :inputs (x observations coefficients)
  :outputs (residuals)
  :gsl-version (1 11)
  :documentation			; FDL
  "Compute the vector of residuals r = y - X c for the observations y,
  coefficients c and matrix of predictor variables X.")

298
299
300
301
;;;;****************************************************************************
;;;; Examples
;;;;****************************************************************************

302
(defun linear-least-squares-univariate-example (&optional (print-steps t))
303
304
  "First example in Section 36.5 of the GSL manual."
  ;; Results not given in manual so not verified yet.
305
306
307
308
309
310
311
312
313
314
315
316
  (let ((x
	 (grid:make-foreign-array
	  'double-float :initial-contents
	  '(1970.0d0 1980.0d0 1990.0d0 2000.0d0)))
	(y
	 (grid:make-foreign-array
	  'double-float :initial-contents
	  '(12.0d0 11.0d0 14.0d0 13.0d0)))
	(w
	 (grid:make-foreign-array
	  'double-float :initial-contents
	  '(0.1d0 0.2d0 0.3d0 0.4d0))))
317
    (multiple-value-bind (c0 c1 cov00 cov01 cov11 chisq)
318
	(linear-fit x y w)
319
320
321
322
323
324
      (when print-steps
	(format t "Best fit: Y = ~8,5f + ~8,5f X~&" c0 c1)
	(format t "Covariance matrix:~&[~12,5f ~12,5f~&~12,5f ~12,5f]~&"
		cov00 cov01 cov01 cov11)
	(format t "Chisq = ~g~&" chisq)
	(loop for i from 0 below (dim0 x)
325
	      do
326
	   (format t "data: ~12,5f ~12,5f ~12,5f~&"
327
328
329
		     (grid:aref x i)
		     (grid:aref y i)
		     (/ (grid:aref w i))))
330
	(loop for i from -30 below 130 by 10 ; don't print everything
331
332
333
334
335
336
	      for
	      xf = (+ (grid:aref x 0)
		      (* (/ i 100)
			 (- (grid:aref x (1- (dim0 x)))
			    (grid:aref x 0))))
	      do
337
338
339
340
341
342
	   (multiple-value-bind (yf yferr)
	       (linear-estimate xf c0 c1 cov00 cov01 cov11)
	     (format t "fit:~6t~g ~g~&" xf yf)
	     (format t "high:~6t~g ~g~&" xf (+ yf yferr))
	     (format t "low:~6t~g ~g~&" xf (- yf yferr)))))
      (values c0 c1 cov00 cov01 cov11 chisq))))
343
344
345
346

(defun mv-linear-least-squares-data ()
  "Generate data for second example in Section 36.5 of the GSL
   manual."
Liam Healy's avatar
Liam Healy committed
347
  (let ((rng (make-random-number-generator +mt19937+ 0)))
348
349
350
351
352
    (loop for x from 1/10 below 2 by 1/10
	  for xd = (coerce x 'double-float)
	  for y0 = (exp xd)
	  for sigma = (* 0.1d0 y0)
	  collect
353
	  (list xd (+ y0 (sample rng :gaussian :sigma sigma)) sigma))))
354

355
(defun linear-least-squares-multivariate-example (data &optional (print-details t))
356
357
358
359
  "Second example in Section 36.5 of the GSL manual.  Returns the
   coefficients of x^0, x^1, x^2 for the best fit, and the chi
   squared."
  (let* ((n (length data))
360
361
362
	 (x (grid:make-foreign-array 'double-float :dimensions (list n 3)))
	 (y (grid:make-foreign-array 'double-float :dimensions n))
	 (w (grid:make-foreign-array 'double-float :dimensions n)))
363
    (loop for i from 0
364
       for row in data do
365
366
367
368
369
       (setf (grid:aref X i 0) 1.0d0
	     (grid:aref X i 1) (first row)
	     (grid:aref X i 2) (expt (first row) 2)
	     (grid:aref y i) (second row)
	     (grid:aref w i) (/ (expt (third row) 2))))
370
    (multiple-value-bind (parameters cov chisq)
371
	(linear-mfit X y 3 w)
372
      (when print-details
373
	(format t "Best fit: Y = ~10,8f + ~10,8f X + ~10,8f X^2~&"
374
		(grid:aref parameters 0) (grid:aref parameters 1) (grid:aref parameters 2))
375
376
	(format t "Covariance matrix:~&")
	(format t "~10,8f ~10,8f ~10,8f~&"
377
		(grid:aref cov 0 0) (grid:aref cov 0 1) (grid:aref cov 0 2))
378
	(format t "~10,8f ~10,8f ~10,8f~&"
379
		(grid:aref cov 1 0) (grid:aref cov 1 1) (grid:aref cov 1 2))
380
	(format t "~10,8f ~10,8f ~10,8f~&"
381
		(grid:aref cov 2 0) (grid:aref cov 2 1) (grid:aref cov 2 2))
382
383
	(format t "Chisq = ~10,6f~&" chisq))
      (values
384
       (grid:aref parameters 0) (grid:aref parameters 1) (grid:aref parameters 2)
385
       chisq))))
386

387
(save-test linear-least-squares
388
389
 (linear-least-squares-univariate-example nil)
 (linear-least-squares-multivariate-example (mv-linear-least-squares-data) nil))