Commit f755beaf authored by Liam Healy's avatar Liam Healy

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.
parent f15ffdda
;; 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")
......
;; 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))
;; 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)))
;; 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))
;; 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))))
;; 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)))))
......@@ -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))))
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment