linear-least-squares.lisp 13.8 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: <2008-12-26 18:38:44EST linear-least-squares.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5 6 7

(in-package :gsl)

8 9
;;; /usr/include/gsl/gsl_fit.h

10 11 12 13 14 15 16 17 18 19
;;;;****************************************************************************
;;;; Linear regression
;;;;****************************************************************************

;;; Error in GSL documentation or in code: vectors must be equal
;;; lengths but strides re not?  Shouldn't (floor length stride) be n?
;;; 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.

liam's avatar
liam committed
20
(defmfun linear-fit (x y &optional (x-stride 1) (y-stride 1))
21
  "gsl_fit_linear"
22 23 24
  (((c-pointer x) :pointer) (x-stride sizet)
   ((c-pointer y) :pointer) (y-stride sizet)
   ((dim0 x) sizet) (c0 :double) (c1 :double)
25 26
   (cov00 :double) (cov01 :double) (cov11 :double)
   (sumsq :double))
27
  :inputs (x y)
28 29 30 31 32 33 34 35 36 37 38 39 40
  :documentation			; FDL
  "Determine the best-fit linear regression coefficients
   and returns as the first two values c0,c1
   of the model Y = c_0 + c_1 X for the dataset
   (x,y), two vectors of equal length with strides
   x-stride and y-stride.  The errors on y are assumed unknown so 
   the variance-covariance matrix for the
   parameters (c0, c1) is estimated from the scatter of the
   points around the best-fit line and returned as the third, fourth
   and fifth values cov00, cov01, cov11.   
   The sum of squares of the residuals from the best-fit line is returned
   as the final value.")

liam's avatar
liam committed
41
(defmfun weighted-linear-fit
42 43
    (x weight y &optional (x-stride 1) (weight-stride 1) (y-stride 1))
  "gsl_fit_wlinear"
44 45 46 47
  (((c-pointer x) :pointer) (x-stride sizet)
   ((c-pointer weight) :pointer) (weight-stride sizet)
   ((c-pointer y) :pointer) (y-stride sizet)
   ((dim0 x) sizet) (c0 :double) (c1 :double)
48 49
   (cov00 :double) (cov01 :double) (cov11 :double)
   (chisq :double))
50
  :inputs (x weight y)
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
  :documentation			; FDL
  "Compute the best-fit linear regression coefficients
   c0, c1 of the model Y = c_0 + c_1 X for the weighted
   dataset (x, y), two vectors of equal length with strides
   x-stride and y-stride, and return as the first two values.
   The vector weight, of the same length
   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
   (cov00, cov01, c0v01) as the next three values.  The weighted sum
   of squares of the residuals from the best-fit line, \chi^2, is
   returned in as the last value.")

liam's avatar
liam committed
66
(defmfun linear-estimate (x c0 c1 cov00 cov01 cov11)
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
  "gsl_fit_linear_est"
  ((x :double) (c0 :double) (c1 :double)
   (cov00 :double) (cov01 :double) (cov11 :double)
   (y :double) (y-error :double))
  :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
;;;;****************************************************************************

liam's avatar
liam committed
82
(defmfun multiplier-fit (x y &optional (x-stride 1) (y-stride 1))
83
  "gsl_fit_mul"
84 85 86
  (((c-pointer x) :pointer) (x-stride sizet)
   ((c-pointer y) :pointer) (y-stride sizet)
   ((dim0 x) sizet) (c1 :double) (cov11 :double)
87
   (sumsq :double))
88
  :inputs (x y)
89 90 91 92 93 94 95 96 97
  :documentation			; FDL
  "The best-fit linear regression coefficient c1 of the model Y = c_1
   X for the datasets (x, y) two vectors of equal length with strides x-stride
   and y-stride.  The errors on y are assumed unknown so the variance of
   the parameter c1 is estimated from the scatter of the points around
   the best-fit line and returned as the the second value.  The sum of
   squares of the residuals from the best-fit line is returned as the
   last value.")

liam's avatar
liam committed
98
(defmfun weighted-multiplier-fit
99 100
    (x weight y &optional (x-stride 1) (weight-stride 1) (y-stride 1))
  "gsl_fit_wmul"
101 102 103 104
  (((c-pointer x) :pointer) (x-stride sizet)
   ((c-pointer weight) :pointer) (weight-stride sizet)
   ((c-pointer y) :pointer) (y-stride sizet)
   ((dim0 x) sizet) (c1 :double) (cov11 :double)
105
   (chisq :double))
106
  :inputs (x weight y)
107 108 109 110 111 112 113 114 115 116 117 118 119
  :documentation			; FDL
  "Compute the best-fit linear regression coefficient
   c1 of the model Y = c_1 X for the weighted datasets
   (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
120
(defmfun multiplier-estimate (x c1 cov11)
121 122 123 124 125 126 127 128 129 130 131 132 133
  "gsl_fit_mul_est"
  ((x :double) (c1 :double) (cov11 :double)
   (y :double) (y-error :double))
  :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
;;;;****************************************************************************

134 135 136 137 138
(defmobject fit-workspace
    "gsl_multifit_linear"
  ((number-of-observations sizet) (number-of-parameters sizet))
  "multi-dimensional root solver with function only"
  "Make a workspace for a multidimensional linear least-squares fit.")
139

liam's avatar
liam committed
140
(defmfun linear-mfit
141 142
    (model observations parameters covariance tolerance workspace)
  "gsl_multifit_linear"
143
  (((mpointer model) :pointer) ((mpointer observations) :pointer)
144
   (tolerance :double)
145
   ((mpointer parameters) :pointer) (covariance :pointer) (chisq :double)
146
   ((mpointer workspace) :pointer))
147 148
  :inputs (model observations)
  :outputs (parameters covariance)
149 150 151 152 153 154
  :documentation			; FDL
  "Compute the best-fit parameters c of the model
   y = X c for the observations y and the matrix of predictor
   variables X.  The variance-covariance matrix of the model
   parameters cov is estimated from the scatter of the observations
   about the best-fit.  The sum of squares of the residuals from the
liam's avatar
liam committed
155
   best-fit, chi^2, is returned.
156 157 158 159 160 161 162

   The best-fit is found by singular value decomposition of the matrix
   X using the preallocated workspace provided. The
   modified Golub-Reinsch SVD algorithm is used, 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.")

liam's avatar
liam committed
163
(defmfun linear-mfit-svd
164 165
    (model observations parameters covariance tolerance workspace)
  "gsl_multifit_linear_svd"
166
  (((mpointer model) :pointer) ((mpointer observations) :pointer)
167
   (tolerance :double)
168 169
   (rank sizet)
   ((mpointer parameters) :pointer) (covariance :pointer) (chisq :double)
170
   ((mpointer workspace) :pointer))
171 172
  :inputs (model observations)
  :outputs (parameters covariance)
liam's avatar
liam committed
173
  :return ((dcref chisq) (scref rank))
174 175 176 177 178
  :documentation			; FDL
  "Compute the best-fit parameters c of the model
   y = X c for the observations y and the matrix of predictor
   variables X.  The variance-covariance matrix of the model
   parameters cov is estimated from the scatter of the observations
179 180
   about the best-fit.  The sum of squares of the residuals
   from the best-fit chi^2, and rank are returned.
181 182 183 184 185 186 187

   The best-fit is found by singular value decomposition of the matrix
   X using the preallocated workspace provided. The
   modified Golub-Reinsch SVD algorithm is used, 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.
   In the this form of the function the components are discarded if the
liam's avatar
liam committed
188 189 190
   ratio of singular values s_i/s_0 falls below the user-specified
   tolerance tolerance, and the effective rank is returned as the
   second value.")
191

liam's avatar
liam committed
192
(defmfun weighted-linear-mfit
193 194
    (model weight observations parameters covariance workspace)
  "gsl_multifit_wlinear"
195 196 197 198 199
  (((mpointer model) :pointer)
   ((mpointer weight) :pointer)
   ((mpointer observations) :pointer)
   ((mpointer parameters) :pointer)
   ((mpointer covariance) :pointer) (chisq :double)
200
   ((mpointer workspace) :pointer))
201 202
  :inputs (model observations)
  :outputs (parameters covariance)
203 204 205 206 207 208
  :documentation			; FDL
  "Compute the best-fit parameters c of the weighted
   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 sum of squares of the residuals from the best-fit,
liam's avatar
liam committed
209
   chi^2, is returned as the last value.
210 211 212 213 214 215

   The best-fit is found by singular value decomposition of the matrix
   model using the preallocated workspace provided. Any
   components which have zero singular value (to machine precision) are
   discarded from the fit.")

liam's avatar
liam committed
216
(defmfun weighted-linear-mfit-svd
217 218
    (model weight observations parameters covariance tolerance workspace)
  "gsl_multifit_wlinear_svd"
219 220 221
  (((mpointer model) :pointer)
   ((mpointer weight) :pointer)
   ((mpointer observations) :pointer)
222
   (tolerance :double)
223 224
   (rank sizet)
   ((mpointer parameters) :pointer) (covariance :pointer) (chisq :double)
225
   ((mpointer workspace) :pointer))
226 227
  :inputs (model weight observations)
  :outputs (parameters covariance)
228
  :return ((dcref chisq) (scref rank))
229 230 231 232
  :documentation			; FDL
  "Compute the best-fit parameters c of the weighted
   model y = X c for the observations y and weights
   and the model matrix X.  The covariance matrix of
233 234
   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
235
   chi^2, is returned as the first value.
236 237 238 239 240 241 242 243 244

   The best-fit is found by singular value decomposition of the matrix
   model using the preallocated workspace provided. 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.")

liam's avatar
liam committed
245
(defmfun multi-linear-estimate (x coefficients covariance)
246
  "gsl_multifit_linear_est"
247 248
  (((mpointer x) :pointer) ((mpointer coefficients) :pointer)
   ((mpointer covariance) :pointer) (y :double) (y-error :double))
249
  :inputs (x coefficients covariance)
250 251 252 253 254 255 256 257 258 259 260 261 262
  :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.")

;;;;****************************************************************************
;;;; Examples
;;;;****************************************************************************

(defun univariate-linear-least-squares-example ()
  "First example in Section 36.5 of the GSL manual."
  ;; Results not given in manual so not verified yet.
263
  (let ((x #m(1970.0d0 1980.0d0 1990.0d0 2000.0d0))
264 265
	 (y #m(12.0d0 11.0d0 14.0d0 13.0d0))
	 (w #m(0.1d0 0.2d0 0.3d0 0.4d0)))
266 267 268 269 270 271 272 273 274
	(multiple-value-bind (c0 c1 cov00 cov01 cov11 chisq)
	    (weighted-linear-fit x w y)
	  (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)
		do
		(format t "~&data: ~12,5f ~12,5f ~12,5f"
275 276 277
			(maref x i)
			(maref y i)
			(/ (maref w i))))
278 279
	  (loop for i from -30 below 130 by 10 ; don't print everything
		for
280
		xf = (+ (maref x 0)
281
			(* (/ i 100)
282 283
			   (- (maref x (1- (dim0 x)))
			      (maref x 0))))
284 285 286 287 288 289
		do
		(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))))
290
	  (fresh-line))))
291 292 293 294

(defun mv-linear-least-squares-data ()
  "Generate data for second example in Section 36.5 of the GSL
   manual."
295
  (let ((rng (make-random-number-generator *mt19937* 0)))
296 297 298 299 300 301
    (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
	  (list xd (+ y0 (gaussian rng sigma)) sigma))))
302 303 304

(defun mv-linear-least-squares-example (data)
  "Second example in Section 36.5 of the GSL manual."
305
  (let* ((n (length data)) chisq
306 307 308 309 310
	 (x (make-marray 'double-float :dimensions (list n 3)))
	 (cov (make-marray 'double-float :dimensions '(3 3)))
	 (y (make-marray 'double-float :dimensions n))
	 (w (make-marray 'double-float :dimensions n))
	 (c (make-marray 'double-float :dimensions 3)))
311
    (loop for i from 0
312 313 314 315 316 317
       for row in data do
       (setf (maref X i 0) 1.0d0
	     (maref X i 1) (first row)
	     (maref X i 2) (expt (first row) 2)
	     (maref y i) (second row)
	     (maref w i) (/ (expt (third row) 2))))
318
    (let ((ws (make-fit-workspace n 3)))
319 320 321
      (setf chisq
	    (weighted-linear-mfit X w y c cov ws)))
    (format t "~&Best fit: Y = ~10,8f + ~10,8f X + ~10,8f X^2"
322
	    (maref c 0) (maref c 1) (maref c 2))
323 324 325
    (format t "~&Covariance matrix:")
    (format
     t "~&~10,8f ~10,8f ~10,8f"
326
     (maref cov 0 0) (maref cov 0 1) (maref cov 0 2))
327 328
    (format
     t "~&~10,8f ~10,8f ~10,8f"
329
     (maref cov 1 0) (maref cov 1 1) (maref cov 1 2))
330 331
    (format
     t "~&~10,8f ~10,8f ~10,8f"
332
     (maref cov 2 0) (maref cov 2 1) (maref cov 2 2))
333
    (format t "~&Chisq = ~10,6f" chisq)))
334 335

;;; (mv-linear-least-squares-example (mv-linear-least-squares-data))