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))