Commit ff9126d9 authored by Liam Healy's avatar Liam Healy
Browse files

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.
parent 8047630a
;; 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))
;; 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))
;; 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))
;; 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)
......
;; 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)
......
;; 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
;; 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)
......
;; 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))))))
......@@ -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+)
......
......@@ -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+)
......
;; 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
......@@ -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+)
......
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