Commit 0d487e0e authored by Liam Healy's avatar Liam Healy

Use GSL tests for random distributions: tdist, laplace

Use GSL tests for random distributions.  Port all tdist tests; some of
laplace.
parent 80387572
;; Definition of GSLL system
;; Liam Healy
;; Time-stamp: <2010-02-07 16:58:10EST gsll.asd>
;; Time-stamp: <2010-04-17 18:47:14EDT gsll.asd>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -160,6 +160,7 @@
((:file "rng-types")
(:file "generators" :depends-on ("rng-types"))
(:file "quasi" :depends-on ("rng-types" "generators"))
(:file "tests")
(:file "gaussian" :depends-on ("rng-types"))
(:file "gaussian-tail" :depends-on ("rng-types"))
(:file "gaussian-bivariate" :depends-on ("rng-types"))
......
;; Mathematical functions
;; Liam Healy, Wed Mar 8 2006 - 22:09
;; Time-stamp: <2009-12-27 09:55:41EST mathematical.lisp>
;; Time-stamp: <2010-04-18 00:22:18EDT mathematical.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -75,7 +75,7 @@
:return ((pmnil cr))
:documentation ; FDL
"Return +1 if x is positive infinity, -1 if negative infinity
nil if finite.")
nil if finite. Some platforms will return only +1 for either sign.")
(defmfun finitep (x)
"gsl_finite" ((x :double))
......
;; Exponential distribution
;; Liam Healy, Sun Sep 17 2006
;; Time-stamp: <2010-01-17 10:25:36EST laplace.lisp>
;; Time-stamp: <2010-04-17 18:49:47EDT laplace.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -73,12 +73,8 @@
;;; Examples and unit test
(save-test laplace
(let ((rng (make-random-number-generator +mt19937+ 0)))
(loop for i from 0 to 10
collect
(sample rng :laplace :a 10.0d0)))
(laplace-pdf 0.0d0 10.0d0)
(laplace-p 1.0d0 2.0d0)
(laplace-q 1.0d0 2.0d0)
(laplace-pinv 0.6967346701436833d0 2.0d0)
(laplace-qinv 0.3032653298563167d0 2.0d0))
(testpdf 'laplace-pdf :laplace :a 2.75d0)
(laplace-p 1.0d0 2.0d0)
(laplace-q 1.0d0 2.0d0)
(laplace-pinv 0.6967346701436833d0 2.0d0)
(laplace-qinv 0.3032653298563167d0 2.0d0))
This diff is collapsed.
;; 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>
;;
;; Copyright 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
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program. If not, see <http://www.gnu.org/licenses/>.
(in-package :gsl)
;; randist/test.c
(defconstant +gslt-BINS+ 100)
(defconstant +gslt-lower-limit+ -5.0d0) ; a
(defconstant +gslt-upper-limit+ 5.0d0) ; b
(defconstant +gslt-bin-size+ ; dx
(/ (- +gslt-upper-limit+ +gslt-lower-limit+) +gslt-BINS+))
(defconstant +initial-number-of-samples+ 100000) ; N
(defun distribution-bin-integral (pdf)
"Compute the CDF for a bin from the PDF."
(let ((cdf (make-array +gslt-BINS+ :element-type 'double-float)))
(dotimes (i +gslt-BINS+ cdf)
(let* ((xtrial (+ +gslt-lower-limit+ (* i +gslt-bin-size+)))
(x (if (< (abs xtrial) 1.0d-10) 0.0d0 xtrial)))
(setf (aref cdf i)
(integration-qags
pdf x (+ x +gslt-bin-size+) 1.0d-16 1.0d-4 1000))))))
(defun bin-samples (count edge number-of-samples distribution &rest distribution-parameters)
"Update the fixnum vectors count and edge with sample random values from the distribution,
and return the mean."
(let ((rng (make-random-number-generator +default-type+))
(total 0.0d0))
(dotimes (i number-of-samples)
(let ((r (apply 'sample rng distribution distribution-parameters)))
(incf total r)
(when (< +gslt-lower-limit+ r +gslt-upper-limit+)
(multiple-value-bind (bin frac)
(floor (- r +gslt-lower-limit+) +gslt-bin-size+)
(if (zerop frac)
(incf (aref edge bin))
(incf (aref count bin)))))))
;; If the bin above is empty, its lower edge hits belong in the lower bin.
(dotimes (i +gslt-BINS+)
(when (and (< (1+ i) +gslt-BINS+) (zerop (aref count i)))
(incf (aref count i) (aref edge (1+ i)))
(setf (aref edge (1+ i)) 0)))
(/ total number-of-samples)))
(defun limits-check (count number-of-samples cdf retry)
"Returns nil if finite and within limits."
(let ((retval nil))
(dotimes (i +gslt-BINS+ retval)
(setf retval
(if (not (finitep (aref cdf i)))
t
(if (not (zerop (aref cdf i)))
(let ((d (- (aref count i) (* number-of-samples (aref cdf i)))))
(and (> (/ d (sqrt (* number-of-samples (aref cdf i)))) 5)
(> d 2)))
(zerop (aref count i)))))
(if (and retry retval) (return retval))
;; Printout out of warning will need
;; (x (+ +gslt-lower-limit+ (* i +gslt-bin-size+)))
(if retval (warn "observed vs. expected")))))
(defun testpdf (pdf distribution &rest distribution-parameters)
"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))))))
(dotimes (i 50)
(let ((mean
(apply 'bin-samples
count edge number-of-samples
distribution distribution-parameters)))
(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,36 +19,17 @@
(in-package :gsl)
(LISP-UNIT:DEFINE-TEST LAPLACE
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST
(LIST 0.005166356198580803d0 -3.942577717493807d0
-8.329510281601332d0 1.1159975704649974d0
-6.2234038148786555d0 -35.04800398421181d0
0.8888158320028845d0 7.161892491969776d0
25.24637780914261d0 7.341651048064451d0
6.54142651602034d0))
(MULTIPLE-VALUE-LIST
(LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
(LOOP FOR I FROM 0 TO 10 COLLECT
(sample rng :laplace :a 10.0d0)))))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 0.05d0)
(MULTIPLE-VALUE-LIST
(LAPLACE-PDF 0.0d0
10.0d0)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 0.6967346701436833d0)
(MULTIPLE-VALUE-LIST (LAPLACE-P 1.0d0 2.0d0)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 0.3032653298563167d0)
(MULTIPLE-VALUE-LIST (LAPLACE-Q 1.0d0 2.0d0)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 1.0d0)
(MULTIPLE-VALUE-LIST
(LAPLACE-PINV
0.6967346701436833d0
2.0d0)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 1.0d0)
(MULTIPLE-VALUE-LIST
(LAPLACE-QINV
0.3032653298563167d0
2.0d0))))
(LISP-UNIT::ASSERT-true (testpdf 'laplace-pdf :laplace :a 2.75d0))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 0.6967346701436833d0)
(MULTIPLE-VALUE-LIST (LAPLACE-P 1.0d0 2.0d0)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 0.3032653298563167d0)
(MULTIPLE-VALUE-LIST (LAPLACE-Q 1.0d0 2.0d0)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 1.0d0)
(MULTIPLE-VALUE-LIST (LAPLACE-PINV 0.6967346701436833d0 2.0d0)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 1.0d0)
(MULTIPLE-VALUE-LIST (LAPLACE-QINV 0.3032653298563167d0 2.0d0))))
This diff is collapsed.
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