From ff9126d9fac76de3851246bfb543d624d3129a38 Mon Sep 17 00:00:00 2001 From: Liam Healy <liam@thinkpad.local> Date: Mon, 24 May 2010 21:11:48 -0400 Subject: [PATCH] New random sample tests; rewrite testpdf Rewrote testpdf so that it takes a function of the random variable, and the list of distribution name and keyword-argument pair. Added new tests levy and exponential-power; several more (beta, chi-squared, exponential) written but commented out because they take an unreasonably long time to complete; these require further investigation. --- random/beta.lisp | 24 ++++++++------- random/chi-squared.lisp | 22 +++++++------- random/exponential-power.lisp | 25 ++++++++++------ random/exponential.lisp | 4 ++- random/laplace.lisp | 4 +-- random/levy.lisp | 32 +++++++++++++------- random/tdist.lisp | 6 ++-- random/tests.lisp | 20 ++++--------- tests/exponential-power.lisp | 36 +++++++++++++---------- tests/laplace.lisp | 3 +- tests/levy.lisp | 55 +++++++++++++++++++---------------- tests/tdist.lisp | 6 ++-- 12 files changed, 132 insertions(+), 105 deletions(-) diff --git a/random/beta.lisp b/random/beta.lisp index 341a52b8..de9f9ab7 100644 --- a/random/beta.lisp +++ b/random/beta.lisp @@ -1,6 +1,6 @@ ;; Beta distribution ;; Liam Healy, Sat Sep 30 2006 -;; Time-stamp: <2010-01-17 10:05:21EST beta.lisp> +;; Time-stamp: <2010-05-24 20:11:32EDT beta.lisp> ;; ;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy ;; Distributed under the terms of the GNU General Public License @@ -68,16 +68,18 @@ Q(x) for the beta distribution with parameters a and b.") ;;; Examples and unit test -(save-test beta - (let ((rng (make-random-number-generator +mt19937+ 0))) - (loop for i from 0 to 10 - collect - (sample rng :beta :a 1.0d0 :b 2.0d0))) - (beta-pdf 0.1d0 1.0d0 2.0d0) - (beta-P 0.1d0 1.0d0 2.0d0) - (beta-Q 0.1d0 1.0d0 2.0d0) - (beta-Pinv 0.19d0 1.0d0 2.0d0) - (beta-Qinv 0.81d0 1.0d0 2.0d0)) +(save-test + beta + ;;(testpdf (lambda (r) (beta-pdf r 2.0d0 3.0d0)) :beta :a 2.0d0 :b 3.0d0) + (let ((rng (make-random-number-generator +mt19937+ 0))) + (loop for i from 0 to 10 + collect + (sample rng :beta :a 1.0d0 :b 2.0d0))) + (beta-pdf 0.1d0 1.0d0 2.0d0) + (beta-P 0.1d0 1.0d0 2.0d0) + (beta-Q 0.1d0 1.0d0 2.0d0) + (beta-Pinv 0.19d0 1.0d0 2.0d0) + (beta-Qinv 0.81d0 1.0d0 2.0d0)) diff --git a/random/chi-squared.lisp b/random/chi-squared.lisp index 17b32a34..fa002b5f 100644 --- a/random/chi-squared.lisp +++ b/random/chi-squared.lisp @@ -1,6 +1,6 @@ ;; Chi-squared distribution ;; Liam Healy, Sat Oct 7 2006 - 16:13 -;; Time-stamp: <2010-01-17 10:07:11EST chi-squared.lisp> +;; Time-stamp: <2010-05-24 20:46:27EDT chi-squared.lisp> ;; ;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy ;; Distributed under the terms of the GNU General Public License @@ -73,15 +73,17 @@ ;;; Examples and unit test (save-test chi-squared - (let ((rng (make-random-number-generator +mt19937+ 0))) - (loop for i from 0 to 10 - collect - (sample rng :chi-squared :nu 10.0d0))) - (chi-squared-pdf 0.5d0 1.0d0) - (chi-squared-P 0.5d0 1.0d0) - (chi-squared-Q 0.5d0 1.0d0) - (chi-squared-Pinv 0.5204998778130463d0 1.0d0) - (chi-squared-Qinv 0.4795001221869537d0 1.0d0)) + (let ((rng (make-random-number-generator +mt19937+ 0))) + (loop for i from 0 to 10 + collect + (sample rng :chi-squared :nu 10.0d0))) + ;; From randist/test.c + ;;(testpdf (lambda (r) (chi-squared-pdf r 13.0d0)) :chi-squared :nu 3.0d0) + (chi-squared-pdf 0.5d0 1.0d0) + (chi-squared-P 0.5d0 1.0d0) + (chi-squared-Q 0.5d0 1.0d0) + (chi-squared-Pinv 0.5204998778130463d0 1.0d0) + (chi-squared-Qinv 0.4795001221869537d0 1.0d0)) diff --git a/random/exponential-power.lisp b/random/exponential-power.lisp index 36f900f1..28ce5969 100644 --- a/random/exponential-power.lisp +++ b/random/exponential-power.lisp @@ -1,6 +1,6 @@ ;; Exponential power distribution ;; Liam Healy, Sat Sep 30 2006 -;; Time-stamp: <2010-01-17 10:11:03EST exponential-power.lisp> +;; Time-stamp: <2010-05-24 20:53:28EDT exponential-power.lisp> ;; ;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy ;; Distributed under the terms of the GNU General Public License @@ -63,11 +63,18 @@ parameters a and b.") ;;; Examples and unit test -(save-test exponential-power - (let ((rng (make-random-number-generator +mt19937+ 0))) - (loop for i from 0 to 10 - collect - (sample rng :exponential-power :a 1.0d0 :b 2.0d0))) - (exponential-power-pdf 0.0d0 1.0d0 2.0d0) - (exponential-power-P 1.0d0 1.0d0 2.0d0) - (exponential-power-Q 1.0d0 1.0d0 2.0d0)) +(save-test + exponential-power + ;; From randist/test.c + (testpdf (lambda (r) (exponential-power-pdf r 3.7d0 0.3d0)) ; exppow0 + :exponential-power :a 3.7d0 :b 0.3d0) + (testpdf (lambda (r) (exponential-power-pdf r 3.7d0 1.0d0)) ; exppow1 + :exponential-power :a 3.7d0 :b 1.0d0) + (testpdf (lambda (r) (exponential-power-pdf r 3.7d0 1.9d0)) ; exppow1a + :exponential-power :a 3.7d0 :b 1.9d0) + (testpdf (lambda (r) (exponential-power-pdf r 3.7d0 2.0d0)) ; exppow2 + :exponential-power :a 3.7d0 :b 2.0d0) + (testpdf (lambda (r) (exponential-power-pdf r 3.7d0 3.5d0)) ; exppow2a + :exponential-power :a 3.7d0 :b 3.5d0) + (testpdf (lambda (r) (exponential-power-pdf r 3.7d0 7.5d0)) ; exppow2b + :exponential-power :a 3.7d0 :b 7.5d0)) diff --git a/random/exponential.lisp b/random/exponential.lisp index cdd84eae..c708d7e1 100644 --- a/random/exponential.lisp +++ b/random/exponential.lisp @@ -1,6 +1,6 @@ ;; Exponential distribution ;; Liam Healy, Sat Sep 2 2006 - 19:04 -;; Time-stamp: <2010-01-17 10:12:00EST exponential.lisp> +;; Time-stamp: <2010-05-24 20:46:09EDT exponential.lisp> ;; ;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy ;; Distributed under the terms of the GNU General Public License @@ -78,6 +78,8 @@ (loop for i from 0 to 10 collect (sample rng :exponential :mu 10.0d0))) + ;; From randist/test.c + ;;(testpdf (lambda (r) (exponential-pdf r 2.0d0)) :exponential :mu 2.0d0) (exponential-pdf 0.0d0 10.0d0) (exponential-p 1.0d0 2.0d0) (exponential-q 1.0d0 2.0d0) diff --git a/random/laplace.lisp b/random/laplace.lisp index 9ee11ea7..84063b84 100644 --- a/random/laplace.lisp +++ b/random/laplace.lisp @@ -1,6 +1,6 @@ ;; Exponential distribution ;; Liam Healy, Sun Sep 17 2006 -;; Time-stamp: <2010-04-17 18:49:47EDT laplace.lisp> +;; Time-stamp: <2010-05-24 19:58:32EDT laplace.lisp> ;; ;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy ;; Distributed under the terms of the GNU General Public License @@ -73,7 +73,7 @@ ;;; Examples and unit test (save-test laplace - (testpdf 'laplace-pdf :laplace :a 2.75d0) + (testpdf (lambda (r) (laplace-pdf r 2.75d0)) :laplace :a 2.75d0) (laplace-p 1.0d0 2.0d0) (laplace-q 1.0d0 2.0d0) (laplace-pinv 0.6967346701436833d0 2.0d0) diff --git a/random/levy.lisp b/random/levy.lisp index 8bf0f646..417bb4ec 100644 --- a/random/levy.lisp +++ b/random/levy.lisp @@ -1,8 +1,8 @@ ;; Levy distribution ;; Liam Healy, Sat Sep 30 2006 -;; Time-stamp: <2010-01-17 10:26:16EST levy.lisp> +;; Time-stamp: <2010-05-24 21:06:37EDT levy.lisp> ;; -;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy +;; Copyright 2006, 2007, 2008, 2009, 2010 Liam M. Healy ;; Distributed under the terms of the GNU General Public License ;; ;; This program is free software: you can redistribute it and/or modify @@ -66,11 +66,23 @@ ;;; Examples and unit test (save-test levy - (let ((rng (make-random-number-generator +mt19937+ 0))) - (loop for i from 0 to 10 - collect - (sample rng :levy :c 1.0d0 :alpha 2.0d0))) - (let ((rng (make-random-number-generator +mt19937+ 0))) - (loop for i from 0 to 10 - collect - (sample rng :levy-skew :c 1.0d0 :alpha 2.0d0 :beta 1.0d0)))) + ;; From randist/test.c + (testpdf (lambda (r) (cauchy-pdf r 5.0d0)) :levy :c 5.0d0 :alpha 1.0d0) ; levy1 + (testpdf (lambda (r) (cauchy-pdf r 5.0d0)) :levy :c 5.0d0 :alpha 1.01d0) ; levy1a + (testpdf (lambda (r) (gaussian-pdf r (* (sqrt 2.0d0) 5.0d0))) + :levy :c 5.0d0 :alpha 2.0d0) ; levy2 + (testpdf (lambda (r) (gaussian-pdf r (* (sqrt 2.0d0) 5.0d0))) + :levy :c 5.0d0 :alpha 1.99d0) ; levy2a + (testpdf (lambda (r) (cauchy-pdf r 5.0d0)) + :levy-skew :c 5.0d0 :alpha 1.0d0 :beta 0.0d0) ; levy_skew1 + (testpdf (lambda (r) (cauchy-pdf r 5.0d0)) + :levy-skew :c 5.0d0 :alpha 1.01d0 :beta 0.0d0) ; levy_skew1a + (testpdf (lambda (r) (gaussian-pdf r (* (sqrt 2.0d0) 5.0d0))) + :levy-skew :c 5.0d0 :alpha 2.0d0 :beta 0.0d0) ; levy_skew2 + (testpdf (lambda (r) (gaussian-pdf r (* (sqrt 2.0d0) 5.0d0))) + :levy-skew :c 5.0d0 :alpha 1.99d0 :beta 0.0d0) ; levy_skew2a + (testpdf (lambda (r) (cauchy-pdf r 5.0d0)) + :levy-skew :c 5.0d0 :alpha 1.01d0 :beta 0.001d0) ; levy_skew1b + (testpdf (lambda (r) (gaussian-pdf r (* (sqrt 2.0d0) 5.0d0))) + :levy-skew :c 5.0d0 :alpha 1.99d0 :beta 0.001d0)) ; levy_skew2b + diff --git a/random/tdist.lisp b/random/tdist.lisp index f49d6a05..16c7ca61 100644 --- a/random/tdist.lisp +++ b/random/tdist.lisp @@ -1,6 +1,6 @@ ;; Tdist distribution ;; Liam Healy, Sat Oct 7 2006 - 16:13 -;; Time-stamp: <2010-04-18 00:27:57EDT tdist.lisp> +;; Time-stamp: <2010-05-24 20:01:35EDT tdist.lisp> ;; ;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy ;; Distributed under the terms of the GNU General Public License @@ -73,8 +73,8 @@ ;;; Examples and unit test (save-test tdist - (testpdf 'tdist-pdf :tdist :nu 1.75d0) - (testpdf 'tdist-pdf :tdist :nu 12.75d0) + (testpdf (lambda (r) (tdist-pdf r 1.75d0)) :tdist :nu 1.75d0) + (testpdf (lambda (r) (tdist-pdf r 12.75d0)) :tdist :nu 12.75d0) ;; From test_tdist in cdf/test.c (tdist-P 0.0d0 1.0d0) diff --git a/random/tests.lisp b/random/tests.lisp index b4589f59..a0a4270f 100644 --- a/random/tests.lisp +++ b/random/tests.lisp @@ -1,6 +1,6 @@ ;; Emulate the GSL tests for random distributions. ;; Liam Healy 2010-04-17 09:44:49EDT tests.lisp -;; Time-stamp: <2010-04-17 18:44:57EDT tests.lisp> +;; Time-stamp: <2010-05-24 20:38:47EDT tests.lisp> ;; ;; Copyright 2010 Liam M. Healy ;; Distributed under the terms of the GNU General Public License @@ -77,27 +77,17 @@ ;; (x (+ +gslt-lower-limit+ (* i +gslt-bin-size+))) (if retval (warn "observed vs. expected"))))) -(defun testpdf (pdf distribution &rest distribution-parameters) +(defun testpdf (pdf-function &rest distribution) "Test the probability density function in the same way that GSL does. If everything is functioning correctly, this will return T." (let ((count (make-array +gslt-BINS+ :element-type 'fixnum :initial-element 0)) (edge (make-array +gslt-BINS+ :element-type 'fixnum :initial-element 0)) (number-of-samples +initial-number-of-samples+) - (cdf (distribution-bin-integral - (lambda (x) - (apply pdf x - ;; It is assumed that the arguments to the - ;; sample function are the same as to the - ;; distribution function except without - ;; keywords. - (loop for (nil value) on distribution-parameters by #'cddr - collect value)))))) + (cdf (distribution-bin-integral pdf-function))) (dotimes (i 50) - (let ((mean - (apply 'bin-samples - count edge number-of-samples - distribution distribution-parameters))) + (let ((mean (apply 'bin-samples count edge number-of-samples distribution))) (unless (finitep mean) (error "mean ~a is not finite" mean)) (if (limits-check count number-of-samples cdf (< (1+ i) 50)) (if (plusp i) (setf number-of-samples (* 2 number-of-samples))) (return t)))))) + diff --git a/tests/exponential-power.lisp b/tests/exponential-power.lisp index 5ffd1668..16616eed 100644 --- a/tests/exponential-power.lisp +++ b/tests/exponential-power.lisp @@ -19,22 +19,26 @@ (in-package :gsl) (LISP-UNIT:DEFINE-TEST EXPONENTIAL-POWER - (LISP-UNIT::ASSERT-NUMERICAL-EQUAL - (LIST - (LIST 0.09469475592777954d0 -0.06229680875327071d0 - 1.183985538537803d0 0.5187626019237904d0 - 0.7053564314063956d0 -0.9033303844569821d0 - -1.6947336289940842d0 -0.4803236108055401d0 - -0.027641736349912214d0 0.6318391856046153d0 - -0.012478875227423025d0)) - (MULTIPLE-VALUE-LIST - (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0))) - (LOOP FOR I FROM 0 TO 10 COLLECT - (sample rng :exponential-power :a 1.0d0 :b 2.0d0))))) - (LISP-UNIT::ASSERT-NUMERICAL-EQUAL - (LIST 0.5641895835477557d0) - (MULTIPLE-VALUE-LIST - (EXPONENTIAL-POWER-PDF 0.0d0 1.0d0 2.0d0))) + ;; From randist/test.c + (lisp-unit::assert-true + (testpdf (lambda (r) (exponential-power-pdf r 3.7d0 0.3d0)) ; exppow0 + :exponential-power :a 3.7d0 :b 0.3d0)) + (lisp-unit::assert-true + (testpdf (lambda (r) (exponential-power-pdf r 3.7d0 1.0d0)) ; exppow1 + :exponential-power :a 3.7d0 :b 1.0d0)) + (lisp-unit::assert-true + (testpdf (lambda (r) (exponential-power-pdf r 3.7d0 1.9d0)) ; exppow1a + :exponential-power :a 3.7d0 :b 1.9d0)) + (lisp-unit::assert-true + (testpdf (lambda (r) (exponential-power-pdf r 3.7d0 2.0d0)) ; exppow2 + :exponential-power :a 3.7d0 :b 2.0d0)) + (lisp-unit::assert-true + (testpdf (lambda (r) (exponential-power-pdf r 3.7d0 3.5d0)) ; exppow2a + :exponential-power :a 3.7d0 :b 3.5d0)) + (lisp-unit::assert-true + (testpdf (lambda (r) (exponential-power-pdf r 3.7d0 7.5d0)) ; exppow2b + :exponential-power :a 3.7d0 :b 7.5d0)) + ;; Automatically converted from cdf/test.c (ASSERT-TO-TOLERANCE (EXPONENTIAL-POWER-P -1000.0d0 0.7d0 1.8d0) 0.0d0 +TEST-TOL6+) (ASSERT-TO-TOLERANCE (EXPONENTIAL-POWER-P -0.1d0 0.7d0 1.8d0) 0.42053490828675155d0 +TEST-TOL0+) (ASSERT-TO-TOLERANCE (EXPONENTIAL-POWER-P -1.d-32 0.7d0 1.8d0) 0.5d0 +TEST-TOL0+) diff --git a/tests/laplace.lisp b/tests/laplace.lisp index 688b55d6..2b6ba5fc 100644 --- a/tests/laplace.lisp +++ b/tests/laplace.lisp @@ -19,7 +19,8 @@ (in-package :gsl) (LISP-UNIT:DEFINE-TEST LAPLACE - (LISP-UNIT::ASSERT-true (testpdf 'laplace-pdf :laplace :a 2.75d0)) + ;; From randist/test.c + (LISP-UNIT::ASSERT-true (testpdf (lambda (r) (laplace-pdf r 2.75d0)) :laplace :a 2.75d0)) ;; Automatically converted from cdf/test_auto.c (ASSERT-TO-TOLERANCE (LAPLACE-P -1.d10 1.3d0) 0.0d0 +TEST-TOL6+) (ASSERT-TO-TOLERANCE (LAPLACE-P -1.d9 1.3d0) 0.0d0 +TEST-TOL6+) diff --git a/tests/levy.lisp b/tests/levy.lisp index a33f9b71..682f1fae 100644 --- a/tests/levy.lisp +++ b/tests/levy.lisp @@ -1,6 +1,6 @@ ;; Regression test LEVY for GSLL, automatically generated ;; -;; Copyright 2009 Liam M. Healy +;; Copyright 2009, 2010 Liam M. Healy ;; Distributed under the terms of the GNU General Public License ;; ;; This program is free software: you can redistribute it and/or modify @@ -19,28 +19,33 @@ (in-package :gsl) (LISP-UNIT:DEFINE-TEST LEVY - (LISP-UNIT::ASSERT-NUMERICAL-EQUAL - (LIST - (LIST 2.6941098332360465d0 -0.29395438644676647d0 - -1.2703401352272083d0 1.0771538640113538d0 - 0.13771218406916103d0 0.9419728438107844d0 - -0.5107134674789159d0 0.1648207853689268d0 - -0.14857899041035147d0 -1.9074885744364487d0 - -2.086195213997167d0)) - (MULTIPLE-VALUE-LIST - (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0))) - (LOOP FOR I FROM 0 TO 10 COLLECT - (sample rng :levy :c 1.0d0 :alpha 2.0d0))))) - (LISP-UNIT::ASSERT-NUMERICAL-EQUAL - (LIST - (LIST 2.6941098332360465d0 -0.2939543864467665d0 - -1.2703401352272083d0 1.0771538640113538d0 - 0.13771218406916097d0 0.9419728438107844d0 - -0.510713467478916d0 0.1648207853689266d0 - -0.14857899041035158d0 -1.907488574436449d0 - -2.086195213997167d0)) - (MULTIPLE-VALUE-LIST - (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0))) - (LOOP FOR I FROM 0 TO 10 COLLECT - (sample rng :levy-skew :c 1.0d0 :alpha 2.0d0 :beta 1.0d0)))))) + ;; From randist/test.c + (lisp-unit::assert-true + (testpdf (lambda (r) (cauchy-pdf r 5.0d0)) :levy :c 5.0d0 :alpha 1.0d0)) ; levy1 + (lisp-unit::assert-true + (testpdf (lambda (r) (cauchy-pdf r 5.0d0)) :levy :c 5.0d0 :alpha 1.01d0)) ; levy1a + (lisp-unit::assert-true + (testpdf (lambda (r) (gaussian-pdf r (* (sqrt 2.0d0) 5.0d0))) + :levy :c 5.0d0 :alpha 2.0d0)) ; levy2 + (lisp-unit::assert-true + (testpdf (lambda (r) (gaussian-pdf r (* (sqrt 2.0d0) 5.0d0))) + :levy :c 5.0d0 :alpha 1.99d0)) ; levy2a + (lisp-unit::assert-true + (testpdf (lambda (r) (cauchy-pdf r 5.0d0)) + :levy-skew :c 5.0d0 :alpha 1.0d0 :beta 0.0d0)) ; levy_skew1 + (lisp-unit::assert-true + (testpdf (lambda (r) (cauchy-pdf r 5.0d0)) + :levy-skew :c 5.0d0 :alpha 1.01d0 :beta 0.0d0)) ; levy_skew1a + (lisp-unit::assert-true + (testpdf (lambda (r) (gaussian-pdf r (* (sqrt 2.0d0) 5.0d0))) + :levy-skew :c 5.0d0 :alpha 2.0d0 :beta 0.0d0)) ; levy_skew2 + (lisp-unit::assert-true + (testpdf (lambda (r) (gaussian-pdf r (* (sqrt 2.0d0) 5.0d0))) + :levy-skew :c 5.0d0 :alpha 1.99d0 :beta 0.0d0)) ; levy_skew2a + (lisp-unit::assert-true + (testpdf (lambda (r) (cauchy-pdf r 5.0d0)) + :levy-skew :c 5.0d0 :alpha 1.01d0 :beta 0.001d0)) ; levy_skew1b + (lisp-unit::assert-true + (testpdf (lambda (r) (gaussian-pdf r (* (sqrt 2.0d0) 5.0d0))) + :levy-skew :c 5.0d0 :alpha 1.99d0 :beta 0.001d0))) ; levy_skew2b diff --git a/tests/tdist.lisp b/tests/tdist.lisp index 2ef7d150..588cc857 100644 --- a/tests/tdist.lisp +++ b/tests/tdist.lisp @@ -19,8 +19,10 @@ (in-package :gsl) (lisp-unit:define-test tdist - (lisp-unit::assert-true (testpdf 'tdist-pdf :tdist :nu 1.75d0)) - (lisp-unit::assert-true (testpdf 'tdist-pdf :tdist :nu 12.75d0)) + ;; From randist/test.c + (lisp-unit::assert-true (testpdf (lambda (r) (tdist-pdf r 1.75d0)) :tdist :nu 1.75d0)) + (lisp-unit::assert-true (testpdf (lambda (r) (tdist-pdf r 12.75d0)) :tdist :nu 12.75d0)) + ;; From cdf/test.c (assert-to-tolerance (tdist-P 0.0d0 1.0d0) 0.5d0 +test-tol6+) (assert-to-tolerance (tdist-P 1d-100 1.0d0) 0.5d0 +test-tol6+) (assert-to-tolerance (tdist-P 0.001d0 1.0d0) 5.00318309780080559d-1 +test-tol6+) -- GitLab