linear-least-squares.lisp 15.2 KB
 liam committed Jan 21, 2008 1 2 ;; Linear least squares, or linear regression ;; Liam Healy <2008-01-21 12:41:46EST linear-least-squares.lisp>  Liam Healy committed Jan 24, 2009 3 ;; Time-stamp: <2009-01-24 09:24:25EST linear-least-squares.lisp>  lhealy committed Jul 24, 2008 4 ;; $Id$  liam committed Jan 21, 2008 5 6 7  (in-package :gsl)  8 9 ;;; /usr/include/gsl/gsl_fit.h  liam committed Jan 21, 2008 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 committed Feb 18, 2008 20 (defmfun linear-fit (x y &optional (x-stride 1) (y-stride 1))  liam committed Jan 21, 2008 21  "gsl_fit_linear"  Liam Healy committed Aug 31, 2008 22 23 24  (((c-pointer x) :pointer) (x-stride sizet) ((c-pointer y) :pointer) (y-stride sizet) ((dim0 x) sizet) (c0 :double) (c1 :double)  liam committed Jan 21, 2008 25 26  (cov00 :double) (cov01 :double) (cov11 :double) (sumsq :double))  Liam Healy committed Dec 07, 2008 27  :inputs (x y)  liam committed Jan 21, 2008 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 committed Feb 18, 2008 41 (defmfun weighted-linear-fit  liam committed Jan 21, 2008 42 43  (x weight y &optional (x-stride 1) (weight-stride 1) (y-stride 1)) "gsl_fit_wlinear"  Liam Healy committed Aug 31, 2008 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)  liam committed Jan 21, 2008 48 49  (cov00 :double) (cov01 :double) (cov11 :double) (chisq :double))  Liam Healy committed Dec 07, 2008 50  :inputs (x weight y)  liam committed Jan 21, 2008 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 committed Feb 18, 2008 66 (defmfun linear-estimate (x c0 c1 cov00 cov01 cov11)  liam committed Jan 21, 2008 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 committed Feb 18, 2008 82 (defmfun multiplier-fit (x y &optional (x-stride 1) (y-stride 1))  liam committed Jan 21, 2008 83  "gsl_fit_mul"  Liam Healy committed Aug 31, 2008 84 85 86  (((c-pointer x) :pointer) (x-stride sizet) ((c-pointer y) :pointer) (y-stride sizet) ((dim0 x) sizet) (c1 :double) (cov11 :double)  liam committed Jan 21, 2008 87  (sumsq :double))  Liam Healy committed Dec 07, 2008 88  :inputs (x y)  liam committed Jan 21, 2008 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 committed Feb 18, 2008 98 (defmfun weighted-multiplier-fit  liam committed Jan 21, 2008 99 100  (x weight y &optional (x-stride 1) (weight-stride 1) (y-stride 1)) "gsl_fit_wmul"  Liam Healy committed Aug 31, 2008 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)  liam committed Jan 21, 2008 105  (chisq :double))  Liam Healy committed Dec 07, 2008 106  :inputs (x weight y)  liam committed Jan 21, 2008 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 committed Feb 18, 2008 120 (defmfun multiplier-estimate (x c1 cov11)  liam committed Jan 21, 2008 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 ;;;;****************************************************************************  Liam Healy committed Dec 25, 2008 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.")  liam committed Jan 21, 2008 139   Liam Healy committed Jan 24, 2009 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 (defun size-array (array-or-size) (if (numberp array-or-size) array-or-size (dim0 array-or-size))) (defun default-covariance (parameters-or-size) (make-marray '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)))  liam committed Feb 18, 2008 155 (defmfun linear-mfit  Liam Healy committed Jan 24, 2009 156 157 158 159 160  (model observations parameters-or-size tolerance &optional (covariance (default-covariance parameters-or-size)) (workspace (default-lls-workspace observations parameters-or-size)) &aux (parameters (vdf parameters-or-size)))  liam committed Jan 21, 2008 161  "gsl_multifit_linear"  Liam Healy committed Aug 31, 2008 162  (((mpointer model) :pointer) ((mpointer observations) :pointer)  liam committed Jan 21, 2008 163  (tolerance :double)  Liam Healy committed Aug 31, 2008 164  ((mpointer parameters) :pointer) (covariance :pointer) (chisq :double)  Liam Healy committed Dec 26, 2008 165  ((mpointer workspace) :pointer))  Liam Healy committed Dec 07, 2008 166 167  :inputs (model observations) :outputs (parameters covariance)  liam committed Jan 21, 2008 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 committed Feb 18, 2008 174  best-fit, chi^2, is returned.  liam committed Jan 21, 2008 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 committed Feb 18, 2008 182 (defmfun linear-mfit-svd  Liam Healy committed Jan 24, 2009 183 184 185 186 187  (model observations parameters-or-size tolerance &optional (covariance (default-covariance parameters-or-size)) (workspace (default-lls-workspace observations parameters-or-size)) &aux (parameters (vdf parameters-or-size)))  liam committed Jan 21, 2008 188  "gsl_multifit_linear_svd"  Liam Healy committed Aug 31, 2008 189  (((mpointer model) :pointer) ((mpointer observations) :pointer)  liam committed Jan 21, 2008 190  (tolerance :double)  Liam Healy committed Aug 31, 2008 191 192  (rank sizet) ((mpointer parameters) :pointer) (covariance :pointer) (chisq :double)  Liam Healy committed Dec 26, 2008 193  ((mpointer workspace) :pointer))  Liam Healy committed Dec 07, 2008 194 195  :inputs (model observations) :outputs (parameters covariance)  liam committed Feb 18, 2008 196  :return ((dcref chisq) (scref rank))  liam committed Jan 21, 2008 197 198 199 200 201  :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  Liam Healy committed Dec 07, 2008 202 203  about the best-fit. The sum of squares of the residuals from the best-fit chi^2, and rank are returned.  liam committed Jan 21, 2008 204 205 206 207 208 209 210  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 committed Feb 18, 2008 211 212 213  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.")  liam committed Jan 21, 2008 214   liam committed Feb 18, 2008 215 (defmfun weighted-linear-mfit  Liam Healy committed Jan 19, 2009 216 217  (model weight observations parameters-or-size &optional  Liam Healy committed Jan 24, 2009 218 219  (covariance (default-covariance parameters-or-size)) (workspace (default-lls-workspace observations parameters-or-size))  Liam Healy committed Jan 19, 2009 220 221  &aux (parameters (vdf parameters-or-size)))  liam committed Jan 21, 2008 222  "gsl_multifit_wlinear"  Liam Healy committed Aug 31, 2008 223 224 225 226 227  (((mpointer model) :pointer) ((mpointer weight) :pointer) ((mpointer observations) :pointer) ((mpointer parameters) :pointer) ((mpointer covariance) :pointer) (chisq :double)  Liam Healy committed Dec 26, 2008 228  ((mpointer workspace) :pointer))  Liam Healy committed Dec 07, 2008 229 230  :inputs (model observations) :outputs (parameters covariance)  Liam Healy committed Jan 19, 2009 231  :return (parameters covariance (dcref chisq))  liam committed Jan 21, 2008 232 233 234 235 236 237  :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 committed Feb 18, 2008 238  chi^2, is returned as the last value.  liam committed Jan 21, 2008 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.")  liam committed Feb 18, 2008 245 (defmfun weighted-linear-mfit-svd  Liam Healy committed Jan 24, 2009 246 247 248 249 250  (model weight observations parameters-or-size tolerance &optional (covariance (default-covariance parameters-or-size)) (workspace (default-lls-workspace observations parameters-or-size)) &aux (parameters (vdf parameters-or-size)))  liam committed Jan 21, 2008 251  "gsl_multifit_wlinear_svd"  Liam Healy committed Aug 31, 2008 252 253 254  (((mpointer model) :pointer) ((mpointer weight) :pointer) ((mpointer observations) :pointer)  liam committed Jan 21, 2008 255  (tolerance :double)  Liam Healy committed Aug 31, 2008 256 257  (rank sizet) ((mpointer parameters) :pointer) (covariance :pointer) (chisq :double)  Liam Healy committed Dec 26, 2008 258  ((mpointer workspace) :pointer))  Liam Healy committed Dec 07, 2008 259 260  :inputs (model weight observations) :outputs (parameters covariance)  liam committed Jan 29, 2008 261  :return ((dcref chisq) (scref rank))  liam committed Jan 21, 2008 262 263 264 265  :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  Liam Healy committed Dec 07, 2008 266 267  the model parameters is computed with the given weights. The weighted sum of squares of the residuals from the best-fit,  liam committed Feb 18, 2008 268  chi^2, is returned as the first value.  liam committed Jan 21, 2008 269 270 271 272 273 274 275 276 277  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 committed Feb 18, 2008 278 (defmfun multi-linear-estimate (x coefficients covariance)  liam committed Jan 21, 2008 279  "gsl_multifit_linear_est"  Liam Healy committed Aug 31, 2008 280 281  (((mpointer x) :pointer) ((mpointer coefficients) :pointer) ((mpointer covariance) :pointer) (y :double) (y-error :double))  Liam Healy committed Dec 07, 2008 282  :inputs (x coefficients covariance)  liam committed Jan 21, 2008 283 284 285 286 287 288 289 290 291 292  :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 ;;;;****************************************************************************  Liam Healy committed Jan 24, 2009 293 (defun linear-least-squares-univariate-example (&optional (print-steps t))  liam committed Jan 21, 2008 294 295  "First example in Section 36.5 of the GSL manual." ;; Results not given in manual so not verified yet.  Liam Healy committed Dec 26, 2008 296  (let ((x #m(1970.0d0 1980.0d0 1990.0d0 2000.0d0))  Liam Healy committed Jan 19, 2009 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324  (y #m(12.0d0 11.0d0 14.0d0 13.0d0)) (w #m(0.1d0 0.2d0 0.3d0 0.4d0))) (multiple-value-bind (c0 c1 cov00 cov01 cov11 chisq) (weighted-linear-fit x w y) (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) do (format t "data: ~12,5f ~12,5f ~12,5f~&" (maref x i) (maref y i) (/ (maref w i)))) (loop for i from -30 below 130 by 10 ; don't print everything for xf = (+ (maref x 0) (* (/ i 100) (- (maref x (1- (dim0 x))) (maref x 0)))) 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))))) (values c0 c1 cov00 cov01 cov11 chisq))))  liam committed Jan 21, 2008 325 326 327 328  (defun mv-linear-least-squares-data () "Generate data for second example in Section 36.5 of the GSL manual."  Liam Healy committed Dec 26, 2008 329  (let ((rng (make-random-number-generator *mt19937* 0)))  liam committed Feb 03, 2008 330 331 332 333 334 335  (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))))  liam committed Jan 21, 2008 336   Liam Healy committed Jan 24, 2009 337 (defun linear-least-squares-multivariate-example (data &optional (print-details t))  Liam Healy committed Jan 19, 2009 338 339 340 341  "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))  Liam Healy committed Dec 26, 2008 342 343  (x (make-marray 'double-float :dimensions (list n 3))) (y (make-marray 'double-float :dimensions n))  Liam Healy committed Jan 19, 2009 344  (w (make-marray 'double-float :dimensions n)))  liam committed Feb 04, 2008 345  (loop for i from 0  Liam Healy committed Nov 30, 2008 346 347 348 349 350 351  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))))  Liam Healy committed Jan 19, 2009 352 353  (multiple-value-bind (parameters cov chisq) (weighted-linear-mfit X w y 3)  Liam Healy committed Jan 24, 2009 354  (when print-details  Liam Healy committed Jan 19, 2009 355 356 357 358 359 360 361 362 363 364 365 366 367  (format t "Best fit: Y = ~10,8f + ~10,8f X + ~10,8f X^2~&" (maref parameters 0) (maref parameters 1) (maref parameters 2)) (format t "Covariance matrix:~&") (format t "~10,8f ~10,8f ~10,8f~&" (maref cov 0 0) (maref cov 0 1) (maref cov 0 2)) (format t "~10,8f ~10,8f ~10,8f~&" (maref cov 1 0) (maref cov 1 1) (maref cov 1 2)) (format t "~10,8f ~10,8f ~10,8f~&" (maref cov 2 0) (maref cov 2 1) (maref cov 2 2)) (format t "Chisq = ~10,6f~&" chisq)) (values (maref parameters 0) (maref parameters 1) (maref parameters 2) chisq))))  liam committed Jan 21, 2008 368   Liam Healy committed Jan 19, 2009 369 (save-test linear-least-squares  Liam Healy committed Jan 24, 2009 370 371  (linear-least-squares-univariate-example nil) (linear-least-squares-multivariate-example (mv-linear-least-squares-data) nil))