From f755beaf84a4ead969036aa4791c42006c6fe3bd Mon Sep 17 00:00:00 2001 From: Liam Healy <liam@thinkpad.local> Date: Mon, 19 Jan 2009 23:34:50 -0500 Subject: [PATCH] Least squares added to tests, make arguments optional Coverted required args to optional args defmfun weighted-linear-mfit and replaced 'parameters with 'parameters-or-size. Added print-steps options for several more examples, and put the newline at the end of the format string. Added tests; CCL 64 bit now fails one of the linear least squares (returns zeros) and the nonlinear least squares (good to 6 or so significant figures). SBCL 64: TOTAL: 1231 assertions passed, 5 failed, 0 execution errors. CCL 64: TOTAL: 1234 assertions passed, 2 failed, 0 execution errors. --- gsll-tests.asd | 4 +- solve-minimize-fit/linear-least-squares.lisp | 126 ++++++++++-------- .../nonlinear-least-squares.lisp | 43 +++--- solve-minimize-fit/roots-one.lisp | 21 +-- tests/linear-least-squares.lisp | 18 +++ tests/nonlinear-least-squares.lisp | 14 ++ tests/roots-one.lisp | 3 + 7 files changed, 146 insertions(+), 83 deletions(-) create mode 100644 tests/linear-least-squares.lisp create mode 100644 tests/nonlinear-least-squares.lisp diff --git a/gsll-tests.asd b/gsll-tests.asd index 2cec6e36..21569e37 100644 --- a/gsll-tests.asd +++ b/gsll-tests.asd @@ -1,6 +1,6 @@ ;; Definition of GSLL system ;; Liam Healy -;; Time-stamp: <2009-01-19 16:40:15EST gsll-tests.asd> +;; Time-stamp: <2009-01-19 22:43:17EST gsll-tests.asd> ;; $Id$ (asdf:defsystem "gsll-tests" @@ -81,6 +81,7 @@ (:file "laplace") (:file "legendre") (:file "levy") + (:file "linear-least-squares") (:file "logarithmic") (:file "logarithm") (:file "logistic") @@ -118,6 +119,7 @@ (:file "monte-carlo") (:file "multinomial") (:file "negative-binomial") + (:file "nonlinear-least-squares") (:file "numerical-differentiation") (:file "numerical-integration") (:file "ode") diff --git a/solve-minimize-fit/linear-least-squares.lisp b/solve-minimize-fit/linear-least-squares.lisp index 2838f1b6..c9881689 100644 --- a/solve-minimize-fit/linear-least-squares.lisp +++ b/solve-minimize-fit/linear-least-squares.lisp @@ -1,6 +1,6 @@ ;; Linear least squares, or linear regression ;; Liam Healy <2008-01-21 12:41:46EST linear-least-squares.lisp> -;; Time-stamp: <2008-12-26 18:38:44EST linear-least-squares.lisp> +;; Time-stamp: <2009-01-19 23:23:11EST linear-least-squares.lisp> ;; $Id$ (in-package :gsl) @@ -189,8 +189,24 @@ tolerance tolerance, and the effective rank is returned as the second value.") +(defun size-array (array-or-size) + (if (numberp array-or-size) + array-or-size + (dim0 array-or-size))) + (defmfun weighted-linear-mfit - (model weight observations parameters covariance workspace) + (model weight observations parameters-or-size + &optional + (covariance + (make-marray + 'double-float + :dimensions + (let ((s (size-array parameters-or-size))) (list s s)))) + (workspace + (make-fit-workspace + (dim0 observations) (size-array parameters-or-size))) + &aux + (parameters (vdf parameters-or-size))) "gsl_multifit_wlinear" (((mpointer model) :pointer) ((mpointer weight) :pointer) @@ -200,6 +216,7 @@ ((mpointer workspace) :pointer)) :inputs (model observations) :outputs (parameters covariance) + :return (parameters covariance (dcref chisq)) :documentation ; FDL "Compute the best-fit parameters c of the weighted model y = X c for the observations y and weights @@ -257,37 +274,38 @@ ;;;; Examples ;;;;**************************************************************************** -(defun univariate-linear-least-squares-example () +(defun univariate-linear-least-squares-example (&optional (print-steps t)) "First example in Section 36.5 of the GSL manual." ;; Results not given in manual so not verified yet. (let ((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))) - (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" - (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)))) - (fresh-line)))) + (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)))) (defun mv-linear-least-squares-data () "Generate data for second example in Section 36.5 of the GSL @@ -300,14 +318,14 @@ collect (list xd (+ y0 (gaussian rng sigma)) sigma)))) -(defun mv-linear-least-squares-example (data) - "Second example in Section 36.5 of the GSL manual." - (let* ((n (length data)) chisq +(defun mv-linear-least-squares-example (data &optional (print-steps t)) + "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)) (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))) + (w (make-marray 'double-float :dimensions n))) (loop for i from 0 for row in data do (setf (maref X i 0) 1.0d0 @@ -315,21 +333,23 @@ (maref X i 2) (expt (first row) 2) (maref y i) (second row) (maref w i) (/ (expt (third row) 2)))) - (let ((ws (make-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" - (maref c 0) (maref c 1) (maref c 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))) + (multiple-value-bind (parameters cov chisq) + (weighted-linear-mfit X w y 3) + (when print-steps + (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)))) -;;; (mv-linear-least-squares-example (mv-linear-least-squares-data)) +(save-test linear-least-squares + (univariate-linear-least-squares-example nil) + (mv-linear-least-squares-example (mv-linear-least-squares-data) nil)) diff --git a/solve-minimize-fit/nonlinear-least-squares.lisp b/solve-minimize-fit/nonlinear-least-squares.lisp index f77a390c..fca32989 100644 --- a/solve-minimize-fit/nonlinear-least-squares.lisp +++ b/solve-minimize-fit/nonlinear-least-squares.lisp @@ -1,6 +1,6 @@ ;; Nonlinear least squares fitting. ;; Liam Healy, 2008-02-09 12:59:16EST nonlinear-least-squares.lisp -;; Time-stamp: <2009-01-03 16:05:48EST nonlinear-least-squares.lisp> +;; Time-stamp: <2009-01-19 22:10:45EST nonlinear-least-squares.lisp> ;; $Id$ (in-package :gsl) @@ -310,11 +310,11 @@ (make-exponent-fit-data :n *number-of-observations* :y - (let ((arr (make-marray 'double-float :dimensions *number-of-observations*))) - (let ((rng (make-random-number-generator *mt19937* 0))) - (dotimes (i *number-of-observations* arr) - (setf (maref arr i) - (+ 1 (* 5 (exp (* -1/10 i))) (gaussian rng 0.1d0)))))) + (let ((arr (make-marray 'double-float :dimensions *number-of-observations*)) + (rng (make-random-number-generator *mt19937* 0))) + (dotimes (i *number-of-observations* arr) + (setf (maref arr i) + (+ 1 (* 5 (exp (* -1/10 i))) (gaussian rng 0.1d0))))) :sigma (make-marray 'double-float :dimensions *number-of-observations* :initial-element 0.1d0)))) @@ -363,7 +363,7 @@ "Find the norm of the fit function f." (euclidean-norm (function-value fit))) -(defun solve-nonlinear-least-squares-example () +(defun solve-nonlinear-least-squares-example (&optional (print-steps t)) (let* ((init #m(1.0d0 0.0d0 0.0d0)) (covariance (make-marray 'double-float @@ -377,9 +377,10 @@ init))) (macrolet ((fitx (i) `(maref (solution fit) ,i)) (err (i) `(sqrt (maref covariance ,i ,i)))) - (format t "~&iter: ~d x = ~15,8f ~15,8f ~15,8f |f(x)|=~7,6g" - 0 (fitx 0) (fitx 1) (fitx 2) - (norm-f fit)) + (when print-steps + (format t "~&iter: ~d x = ~15,8f ~15,8f ~15,8f |f(x)|=~7,6g" + 0 (fitx 0) (fitx 1) (fitx 2) + (norm-f fit))) (loop for iter from 0 below 25 until (and (plusp iter) @@ -387,20 +388,20 @@ do (iterate fit) (ls-covariance (jacobian fit) 0.0d0 covariance) - (format t "~&iter: ~d x = ~15,8f ~15,8f ~15,8f |f(x)|=~7,6g" - (1+ iter) (fitx 0) (fitx 1) (fitx 2) - (norm-f fit)) + (when print-steps + (format t "~&iter: ~d x = ~15,8f ~15,8f ~15,8f |f(x)|=~7,6g" + (1+ iter) (fitx 0) (fitx 1) (fitx 2) + (norm-f fit))) finally (let* ((chi (norm-f fit)) (dof (- *number-of-observations* *number-of-parameters*)) (c (max 1.0d0 (/ chi (sqrt dof))))) - (format t "~&chisq/dof = ~g" (/ (expt chi 2) dof)) - (format t "~&A = ~,5f +/- ~,5f" (fitx 0) (* c (err 0))) - (format t "~&lambda = ~,5f +/- ~,5f" (fitx 1) (* c (err 1))) - (format t "~&b = ~,5f +/- ~,5f" (fitx 2) (* c (err 2))) + (when print-steps + (format t "~&chisq/dof = ~g" (/ (expt chi 2) dof)) + (format t "~&A = ~,5f +/- ~,5f" (fitx 0) (* c (err 0))) + (format t "~&lambda = ~,5f +/- ~,5f" (fitx 1) (* c (err 1))) + (format t "~&b = ~,5f +/- ~,5f" (fitx 2) (* c (err 2)))) (return (list (fitx 0) (fitx 1) (fitx 2)))))))) -;;; Run example: -;;; (nlls-setup) -;;; (solve-nonlinear-least-squares-example) -;;; (5.045357801443204d0 0.10404905892045835d0 1.0192487061031013d0) +(save-test nonlinear-least-squares + (progn (nlls-setup) (solve-nonlinear-least-squares-example nil))) diff --git a/solve-minimize-fit/roots-one.lisp b/solve-minimize-fit/roots-one.lisp index 0172c0f3..f65c1d1e 100644 --- a/solve-minimize-fit/roots-one.lisp +++ b/solve-minimize-fit/roots-one.lisp @@ -1,6 +1,6 @@ ;; One-dimensional root solver. ;; Liam Healy -;; Time-stamp: <2009-01-19 16:31:35EST roots-one.lisp> +;; Time-stamp: <2009-01-19 22:04:19EST roots-one.lisp> ;; $Id$ (in-package :gsl) @@ -330,13 +330,14 @@ (def-single-function quadratic) -(defun roots-one-example () +(defun roots-one-example (&optional (print-steps t)) "Solving a quadratic, the example given in Sec. 32.10 of the GSL manual." (let ((max-iter 50) (solver (make-one-dimensional-root-solver-f *brent-fsolver* quadratic 0.0d0 5.0d0))) - (format t "~&iter ~6t [lower ~24tupper] ~36troot ~44terr ~54terr(est)") + (when print-steps + (format t "iter ~6t [lower ~24tupper] ~36troot ~44terr ~54terr(est)~&")) (loop for iter from 0 for root = (solution solver) for lower = (fsolver-lower solver) @@ -345,10 +346,12 @@ while (and (< iter max-iter) (not (root-test-interval lower upper 0.0d0 0.001d0))) do - (format t "~&~d~6t~10,6f~18t~10,6f~28t~12,9f ~44t~10,4g ~10,4g" - iter lower upper - root (- root (sqrt 5.0d0)) - (- upper lower))))) + (when print-steps + (format t "~d~6t~10,6f~18t~10,6f~28t~12,9f ~44t~10,4g ~10,4g~&" + iter lower upper + root (- root (sqrt 5.0d0)) + (- upper lower))) + finally (return root)))) ;;; Because def-solver-functions and def-single-function bind a symbol ;;; of the same name as the first function, and we want both to run, @@ -377,4 +380,6 @@ iter root (- root (sqrt 5.0d0)) (- root oldroot))) finally (return root)))) -(save-test roots-one (roots-one-fdf-example nil)) +(save-test roots-one + (roots-one-example nil) + (roots-one-fdf-example nil)) diff --git a/tests/linear-least-squares.lisp b/tests/linear-least-squares.lisp new file mode 100644 index 00000000..e6b5f226 --- /dev/null +++ b/tests/linear-least-squares.lisp @@ -0,0 +1,18 @@ +;; Regression test LINEAR-LEAST-SQUARES for GSLL, automatically generated + +(in-package :gsl) + +(LISP-UNIT:DEFINE-TEST LINEAR-LEAST-SQUARES + (LISP-UNIT::ASSERT-NUMERICAL-EQUAL + (LIST -106.59999999999998d0 0.05999999999999999d0 + 39601.99999999999d0 -19.9d0 + 0.009999999999999998d0 0.8d0) + (MULTIPLE-VALUE-LIST + (UNIVARIATE-LINEAR-LEAST-SQUARES-EXAMPLE NIL))) + (LISP-UNIT::ASSERT-NUMERICAL-EQUAL + (LIST 1.1824632487186013d0 0.1845715900137661d0 + 1.3031038153096723d0 16.006137036426168d0) + (MULTIPLE-VALUE-LIST + (MV-LINEAR-LEAST-SQUARES-EXAMPLE + (MV-LINEAR-LEAST-SQUARES-DATA) NIL)))) + diff --git a/tests/nonlinear-least-squares.lisp b/tests/nonlinear-least-squares.lisp new file mode 100644 index 00000000..0dcac436 --- /dev/null +++ b/tests/nonlinear-least-squares.lisp @@ -0,0 +1,14 @@ +;; Regression test NONLINEAR-LEAST-SQUARES for GSLL, automatically generated + +(in-package :gsl) + +(LISP-UNIT:DEFINE-TEST NONLINEAR-LEAST-SQUARES + (LISP-UNIT::ASSERT-NUMERICAL-EQUAL + (LIST + (LIST 5.045357801443204d0 0.10404905892045835d0 + 1.0192487061031013d0)) + (MULTIPLE-VALUE-LIST + (PROGN + (NLLS-SETUP) + (SOLVE-NONLINEAR-LEAST-SQUARES-EXAMPLE NIL))))) + diff --git a/tests/roots-one.lisp b/tests/roots-one.lisp index 8cd09628..961282de 100644 --- a/tests/roots-one.lisp +++ b/tests/roots-one.lisp @@ -3,6 +3,9 @@ (in-package :gsl) (LISP-UNIT:DEFINE-TEST ROOTS-ONE + (LISP-UNIT::ASSERT-NUMERICAL-EQUAL + (LIST 2.2360634081902244d0) + (MULTIPLE-VALUE-LIST (ROOTS-ONE-EXAMPLE NIL))) (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 2.236067977499978d0) (MULTIPLE-VALUE-LIST (ROOTS-ONE-FDF-EXAMPLE NIL)))) -- GitLab