linear-least-squares.lisp 18 KB
Newer Older
1
2
;; Linear least squares, or linear regression
;; Liam Healy <2008-01-21 12:41:46EST linear-least-squares.lisp>
3
;; Time-stamp: <2018-07-14 22:42:07EDT linear-least-squares.lisp>
4
;;
5
;; Copyright 2008, 2009, 2010, 2011, 2012, 2013, 2016, 2018 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
  :documentation			; GSL texi partial
55
  "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
   (cov00, cov01, c0v01) as the next three values.  The weighted or
   unweighted sum of squares of the residuals from the best-fit line,
Liam M. Healy's avatar
Liam M. Healy committed
67
68
   \chi^2, is returned as the last value.

69
70
   Returns: c0, c1, cov00, cov01, cov11, sumsq.
   Returns: intercept, slope, intercept variance, covariance, slope variance, sum square of residuals.")
71

liam's avatar
liam committed
72
(defmfun linear-estimate (x c0 c1 cov00 cov01 cov11)
73
74
75
  "gsl_fit_linear_est"
  ((x :double) (c0 :double) (c1 :double)
   (cov00 :double) (cov01 :double) (cov11 :double)
76
   (y (:pointer :double)) (y-error (:pointer :double)))
77
  :documentation			; GSL texi partial
78
79
80
81
82
83
84
85
86
87
  "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
;;;;****************************************************************************

88
89
90
(defmfun multiplier-fit
    (x y &optional weight (x-stride 1) (y-stride 1) (weight-stride 1))
  ("gsl_fit_mul" "gsl_fit_wmul")
91
92
93
  ((((grid:foreign-pointer x) :pointer) (x-stride :sizet)
    ((grid:foreign-pointer y) :pointer) (y-stride :sizet)
    ((dim0 x) :sizet)
94
95
    (c1 (:pointer :double)) (cov11 (:pointer :double))
    (sumsq (:pointer :double)))
96
97
98
99
   (((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)
100
101
    (c1 (:pointer :double)) (cov11 (:pointer :double))
    (sumsq (:pointer :double))))
102
103
  :inputs (x y weight)
  :switch (weight weight-stride)
104
  :documentation			; GSL texi
105
  "Compute the best-fit linear regression coefficient
106
   c1 of the model Y = c_1 X for the weighted or unweighted datasets
107
108
109
110
111
112
113
114
115
116
   (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
117
(defmfun multiplier-estimate (x c1 cov11)
118
119
  "gsl_fit_mul_est"
  ((x :double) (c1 :double) (cov11 :double)
120
   (y (:pointer :double)) (y-error (:pointer :double)))
121
  :documentation			; GSL texi
122
123
124
125
126
127
128
129
130
  "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
;;;;****************************************************************************

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

138
139
140
141
142
143
(defun size-array (array-or-size)
  (if (numberp array-or-size)
      array-or-size
      (dim0 array-or-size)))

(defun default-covariance (parameters-or-size)
144
  (grid:make-foreign-array
145
146
147
148
149
150
151
152
   '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)))

153
154
155
(export 'linear-mfit)
(defun linear-mfit
    (model observations parameters-or-size
156
     &optional weight tolerance
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
     (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
184
    (model observations parameters-or-size
185
	   &optional
186
	   weight
187
188
	   (covariance (default-covariance parameters-or-size))
	   (workspace (default-lls-workspace observations parameters-or-size))
189
190
	   &aux
	   (parameters (vdf parameters-or-size)))
191
  ("gsl_multifit_linear" "gsl_multifit_wlinear")
192
193
194
195
  ((((mpointer model) :pointer)
    ((mpointer observations) :pointer)
    ((mpointer parameters) :pointer)
    ((mpointer covariance) :pointer)
196
    (chisq (:pointer :double))
197
198
199
200
201
    ((mpointer workspace) :pointer))
   (((mpointer model) :pointer)
    ((mpointer weight) :pointer)
    ((mpointer observations) :pointer)
    ((mpointer parameters) :pointer)
202
    ((mpointer covariance) :pointer)
203
    (chisq (:pointer :double))
204
    ((mpointer workspace) :pointer)))
205
  :inputs (model weight observations)
206
  :outputs (parameters covariance)
207
  :switch (weight)
208
  :return (parameters covariance (cffi:mem-ref chisq :double))
209
210
  :export nil
  :index linear-mfit
211
  :documentation			; GSL texi
212
213
  "Compute the best-fit parameters c of the weighted or unweighted
   model y = X c for the observations y and optional weights
214
215
216
   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
217
   chi^2, is returned as the last value.
218
219

   The best-fit is found by singular value decomposition of the matrix
220
221
222
223
224
   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.")
225

226
#-gsl2 ; GSL v1.x
227
228
229
(defmfun linear-mfit-svd
    (model observations parameters-or-size tolerance
	   &optional weight
230
231
232
	   (covariance (default-covariance parameters-or-size))
	   (workspace (default-lls-workspace observations parameters-or-size))
	   &aux (parameters (vdf parameters-or-size)))
233
  ("gsl_multifit_linear_svd" "gsl_multifit_wlinear_svd")
234
235
  ((((mpointer model) :pointer)
    ((mpointer observations) :pointer)
236
    (tolerance :double)
237
    (rank (:pointer :sizet))
238
239
    ((mpointer parameters) :pointer)
    ((mpointer covariance) :pointer)
240
    (chisq (:pointer :double))
241
242
243
244
245
    ((mpointer workspace) :pointer))
   (((mpointer model) :pointer)
    ((mpointer weight) :pointer)
    ((mpointer observations) :pointer)
    (tolerance :double)
246
    (rank (:pointer :sizet))
247
248
    ((mpointer parameters) :pointer)
    ((mpointer covariance) :pointer)
249
    (chisq (:pointer :double))
250
    ((mpointer workspace) :pointer)))
251
252
  :inputs (model weight observations)
  :outputs (parameters covariance)
253
  :switch (weight)
254
  :return ((cffi:mem-ref chisq :double) (cffi:mem-ref rank ':sizet))
255
256
  :export nil
  :index linear-mfit
257
  :documentation			; GSL texi
Liam M. Healy's avatar
Liam M. Healy committed
258
  "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.
259

Liam M. Healy's avatar
Liam M. Healy committed
260
261
   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, 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.")

262
#+gsl2 ; GSL v2.x
Liam M. Healy's avatar
Liam M. Healy committed
263
264
265
266
267
268
(defmfun linear-mfit-svd
    (model observations parameters-or-size tolerance
	   &optional weight
	   (covariance (default-covariance parameters-or-size))
	   (workspace (default-lls-workspace observations parameters-or-size))
	   &aux (parameters (vdf parameters-or-size)))
269
  ("gsl_multifit_linear_tsvd" "gsl_multifit_wlinear_tsvd")
Liam M. Healy's avatar
Liam M. Healy committed
270
  ((((mpointer model) :pointer)
271
272
273
274
275
276
    ((mpointer observations) :pointer)
    (tolerance :double)
    (rank (:pointer :sizet))
    ((mpointer parameters) :pointer)
    ((mpointer covariance) :pointer)
    (chisq (:pointer :double))
Liam M. Healy's avatar
Liam M. Healy committed
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
    ((mpointer workspace) :pointer))
   (((mpointer model) :pointer)
    ((mpointer weight) :pointer)
    ((mpointer observations) :pointer)
    (tolerance :double)
    (rank (:pointer :sizet))
    ((mpointer parameters) :pointer)
    ((mpointer covariance) :pointer)
    (chisq (:pointer :double))
    ((mpointer workspace) :pointer)))
  :inputs (model weight observations)
  :outputs (parameters covariance)
  :switch (weight)
  :return ((cffi:mem-ref chisq :double) (cffi:mem-ref rank ':sizet))
  :export nil
  :index linear-mfit
  :documentation			; GSL texi
  "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.

   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, 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.")
297

liam's avatar
liam committed
298
(defmfun multi-linear-estimate (x coefficients covariance)
299
  "gsl_multifit_linear_est"
300
  (((mpointer x) :pointer) ((mpointer coefficients) :pointer)
301
302
   ((mpointer covariance) :pointer)
   (y (:pointer :double)) (y-error (:pointer :double)))
303
  :inputs (x coefficients covariance)
304
  :documentation			; GSL texi
305
306
307
308
309
  "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.")

310
311
312
313
(defmfun multi-linear-residuals
    (x observations coefficients
       &optional
       (residuals
314
	(grid:make-foreign-array 'double-float :dimensions (grid:dimensions observations))))
315
316
317
318
319
320
  "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)
321
  :documentation			; GSL texi
322
323
324
  "Compute the vector of residuals r = y - X c for the observations y,
  coefficients c and matrix of predictor variables X.")

325
326
327
328
;;;;****************************************************************************
;;;; Examples
;;;;****************************************************************************

329
(defun linear-least-squares-univariate-example (&optional (print-steps t))
330
331
  "First example in Section 36.5 of the GSL manual."
  ;; Results not given in manual so not verified yet.
332
333
334
335
336
337
338
339
340
341
342
343
  (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))))
344
    (multiple-value-bind (c0 c1 cov00 cov01 cov11 chisq)
345
	(linear-fit x y w)
346
347
348
349
350
351
      (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)
352
	      do
353
	   (format t "data: ~12,5f ~12,5f ~12,5f~&"
354
355
356
		     (grid:aref x i)
		     (grid:aref y i)
		     (/ (grid:aref w i))))
357
	(loop for i from -30 below 130 by 10 ; don't print everything
358
359
360
361
362
363
	      for
	      xf = (+ (grid:aref x 0)
		      (* (/ i 100)
			 (- (grid:aref x (1- (dim0 x)))
			    (grid:aref x 0))))
	      do
364
365
366
367
368
369
	   (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))))
370
371
372
373

(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
374
  (let ((rng (make-random-number-generator +mt19937+ 0)))
375
376
377
378
379
    (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
380
	  (list xd (+ y0 (sample rng :gaussian :sigma sigma)) sigma))))
381

382
(defun linear-least-squares-multivariate-example (data &optional (print-details t))
383
384
385
386
  "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))
387
388
389
	 (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)))
390
    (loop for i from 0
391
       for row in data do
392
393
394
395
396
       (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))))
397
    (multiple-value-bind (parameters cov chisq)
398
	(linear-mfit X y 3 w)
399
      (when print-details
400
	(format t "Best fit: Y = ~10,8f + ~10,8f X + ~10,8f X^2~&"
401
		(grid:aref parameters 0) (grid:aref parameters 1) (grid:aref parameters 2))
402
403
	(format t "Covariance matrix:~&")
	(format t "~10,8f ~10,8f ~10,8f~&"
404
		(grid:aref cov 0 0) (grid:aref cov 0 1) (grid:aref cov 0 2))
405
	(format t "~10,8f ~10,8f ~10,8f~&"
406
		(grid:aref cov 1 0) (grid:aref cov 1 1) (grid:aref cov 1 2))
407
	(format t "~10,8f ~10,8f ~10,8f~&"
408
		(grid:aref cov 2 0) (grid:aref cov 2 1) (grid:aref cov 2 2))
409
410
	(format t "Chisq = ~10,6f~&" chisq))
      (values
411
       (grid:aref parameters 0) (grid:aref parameters 1) (grid:aref parameters 2)
412
       chisq))))
413

414
(save-test linear-least-squares
415
416
 (linear-least-squares-univariate-example nil)
 (linear-least-squares-multivariate-example (mv-linear-least-squares-data) nil))