linear-least-squares.lisp 14.3 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-25 10:35:32EST linear-least-squares.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5 6 7 8 9 10 11 12 13 14 15 16 17

(in-package :gsl)

;;;;****************************************************************************
;;;; 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
18
(defmfun linear-fit (x y &optional (x-stride 1) (y-stride 1))
19
  "gsl_fit_linear"
20 21 22
  (((c-pointer x) :pointer) (x-stride sizet)
   ((c-pointer y) :pointer) (y-stride sizet)
   ((dim0 x) sizet) (c0 :double) (c1 :double)
23 24
   (cov00 :double) (cov01 :double) (cov11 :double)
   (sumsq :double))
25
  :inputs (x y)
26 27 28 29 30 31 32 33 34 35 36 37 38
  :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
39
(defmfun weighted-linear-fit
40 41
    (x weight y &optional (x-stride 1) (weight-stride 1) (y-stride 1))
  "gsl_fit_wlinear"
42 43 44 45
  (((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)
46 47
   (cov00 :double) (cov01 :double) (cov11 :double)
   (chisq :double))
48
  :inputs (x weight y)
49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
  :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
64
(defmfun linear-estimate (x c0 c1 cov00 cov01 cov11)
65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
  "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
80
(defmfun multiplier-fit (x y &optional (x-stride 1) (y-stride 1))
81
  "gsl_fit_mul"
82 83 84
  (((c-pointer x) :pointer) (x-stride sizet)
   ((c-pointer y) :pointer) (y-stride sizet)
   ((dim0 x) sizet) (c1 :double) (cov11 :double)
85
   (sumsq :double))
86
  :inputs (x y)
87 88 89 90 91 92 93 94 95
  :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
96
(defmfun weighted-multiplier-fit
97 98
    (x weight y &optional (x-stride 1) (weight-stride 1) (y-stride 1))
  "gsl_fit_wmul"
99 100 101 102
  (((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)
103
   (chisq :double))
104
  :inputs (x weight y)
105 106 107 108 109 110 111 112 113 114 115 116 117
  :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
118
(defmfun multiplier-estimate (x c1 cov11)
119 120 121 122 123 124 125 126 127 128 129 130 131
  "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
;;;;****************************************************************************

132 133 134 135 136 137 138 139
#|
(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.")
|#

140
(defgo-s (fit-workspace number-of-observations number-of-parameters)
141
	 allocate-fit-workspace free-fit-workspace nil 2)
142

liam's avatar
liam committed
143
(defmfun allocate-fit-workspace (number-of-observations number-of-parameters)
144
  "gsl_multifit_linear_alloc"
145
  ((number-of-observations sizet) (number-of-parameters sizet))
146
  :c-return :pointer
147
  :index (letm fit-workspace)
148 149 150
  :documentation			; FDL
  "Allocate a workspace for fitting a linear model.")

liam's avatar
liam committed
151
(defmfun free-fit-workspace (pointer)
152 153 154
  "gsl_multifit_linear_free"
  ((pointer :pointer))
  :c-return :void
155
  :index (letm fit-workspace)
156 157 158
  :documentation			; FDL
  "Free the memory associate with the workspace.")

liam's avatar
liam committed
159
(defmfun linear-mfit
160 161
    (model observations parameters covariance tolerance workspace)
  "gsl_multifit_linear"
162
  (((mpointer model) :pointer) ((mpointer observations) :pointer)
163
   (tolerance :double)
164
   ((mpointer parameters) :pointer) (covariance :pointer) (chisq :double)
165
   (workspace :pointer))
166 167
  :inputs (model observations)
  :outputs (parameters covariance)
168 169 170 171 172 173
  :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
174
   best-fit, chi^2, is returned.
175 176 177 178 179 180 181

   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
182
(defmfun linear-mfit-svd
183 184
    (model observations parameters covariance tolerance workspace)
  "gsl_multifit_linear_svd"
185
  (((mpointer model) :pointer) ((mpointer observations) :pointer)
186
   (tolerance :double)
187 188
   (rank sizet)
   ((mpointer parameters) :pointer) (covariance :pointer) (chisq :double)
189
   (workspace :pointer))
190 191
  :inputs (model observations)
  :outputs (parameters covariance)
liam's avatar
liam committed
192
  :return ((dcref chisq) (scref rank))
193 194 195 196 197
  :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
198 199
   about the best-fit.  The sum of squares of the residuals
   from the best-fit chi^2, and rank are returned.
200 201 202 203 204 205 206

   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
207 208 209
   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.")
210

liam's avatar
liam committed
211
(defmfun weighted-linear-mfit
212 213
    (model weight observations parameters covariance workspace)
  "gsl_multifit_wlinear"
214 215 216 217 218
  (((mpointer model) :pointer)
   ((mpointer weight) :pointer)
   ((mpointer observations) :pointer)
   ((mpointer parameters) :pointer)
   ((mpointer covariance) :pointer) (chisq :double)
219
   (workspace :pointer))
220 221
  :inputs (model observations)
  :outputs (parameters covariance)
222 223 224 225 226 227
  :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
228
   chi^2, is returned as the last value.
229 230 231 232 233 234

   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
235
(defmfun weighted-linear-mfit-svd
236 237
    (model weight observations parameters covariance tolerance workspace)
  "gsl_multifit_wlinear_svd"
238 239 240
  (((mpointer model) :pointer)
   ((mpointer weight) :pointer)
   ((mpointer observations) :pointer)
241
   (tolerance :double)
242 243
   (rank sizet)
   ((mpointer parameters) :pointer) (covariance :pointer) (chisq :double)
244
   (workspace :pointer))
245 246
  :inputs (model weight observations)
  :outputs (parameters covariance)
247
  :return ((dcref chisq) (scref rank))
248 249 250 251
  :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
252 253
   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
254
   chi^2, is returned as the first value.
255 256 257 258 259 260 261 262 263

   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
264
(defmfun multi-linear-estimate (x coefficients covariance)
265
  "gsl_multifit_linear_est"
266 267
  (((mpointer x) :pointer) ((mpointer coefficients) :pointer)
   ((mpointer covariance) :pointer) (y :double) (y-error :double))
268
  :inputs (x coefficients covariance)
269 270 271 272 273 274 275 276 277 278 279 280 281
  :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.
282 283 284
  (letm ((x #m(1970.0d0 1980.0d0 1990.0d0 2000.0d0))
	 (y #m(12.0d0 11.0d0 14.0d0 13.0d0))
	 (w #m(0.1d0 0.2d0 0.3d0 0.4d0)))
285 286 287 288 289 290 291 292 293
	(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"
294 295 296
			(maref x i)
			(maref y i)
			(/ (maref w i))))
297 298
	  (loop for i from -30 below 130 by 10 ; don't print everything
		for
299
		xf = (+ (maref x 0)
300
			(* (/ i 100)
301 302
			   (- (maref x (1- (dim0 x)))
			      (maref x 0))))
303 304 305 306 307 308
		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))))
309
	  (fresh-line))))
310 311 312 313

(defun mv-linear-least-squares-data ()
  "Generate data for second example in Section 36.5 of the GSL
   manual."
314 315 316 317 318 319 320
  (letm ((rng (random-number-generator *mt19937* 0)))
    (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))))
321 322 323

(defun mv-linear-least-squares-example (data)
  "Second example in Section 36.5 of the GSL manual."
324
  (letm ((n (length data)) chisq
325 326 327 328 329
	 (x (make-array* 'double-float :dimensions (list n 3)))
	 (cov (make-array* 'double-float :dimensions '(3 3)))
	 (y (make-array* 'double-float :dimensions n))
	 (w (make-array* 'double-float :dimensions n))
	 (c (make-array* 'double-float :dimensions 3)))
330
    (loop for i from 0
331 332 333 334 335 336
       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))))
337 338 339 340
    (letm ((ws (fit-workspace n 3)))
      (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"
341
	    (maref c 0) (maref c 1) (maref c 2))
342 343 344
    (format t "~&Covariance matrix:")
    (format
     t "~&~10,8f ~10,8f ~10,8f"
345
     (maref cov 0 0) (maref cov 0 1) (maref cov 0 2))
346 347
    (format
     t "~&~10,8f ~10,8f ~10,8f"
348
     (maref cov 1 0) (maref cov 1 1) (maref cov 1 2))
349 350
    (format
     t "~&~10,8f ~10,8f ~10,8f"
351
     (maref cov 2 0) (maref cov 2 1) (maref cov 2 2))
352
    (format t "~&Chisq = ~10,6f" chisq)))
353 354

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