Commit d0f60ec2 authored by Mirko Vukovic's avatar Mirko Vukovic
Browse files

Merge branch 'x16-rng-chisq-test' into mv

parents 907f96b4 fcb745cf
;;; Code for calculating the chi^2 of a sampled distribution It is
;;; modeled after GSLL's code in test.lisp in the same directory, but
;;; slightly modified. I am not super happy with the aesthetics in
;;; GSLL's code (not very functional), but is is good enough. The
;;; code here is slightly improved, but still not where it ought to
;;; be.
(in-package :gsll)
(defun distribution-bin-integral-1 (pdf bin-count)
"Compute the CDF for a bin from the PDF
Return the CDF array and its sum which can be used to normalize it
Slightly modified from `distribution-bin-integral' in test.lisp"
(let ((cdf (make-array bin-count
:element-type 'double-float
:initial-element 0d0))
(total 0d0))
(dotimes (i bin-count (values cdf total))
(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))
(incf total (aref cdf i))))))
(defun bin-samples-1 (number-of-bins number-of-samples distribution-identifier
&rest distribution-parameters)
"Return the count in each bin and the number of samples in the range
between lower and upper limits The latter can be used for normalizing
purposes
Slightly modified from `distribution-bin-integral' in test.lisp"
(let ((rng (make-random-number-generator +default-type+))
(in-range 0)
(count (make-array number-of-bins
:element-type 'fixnum
:initial-element 0)))
(dotimes (i number-of-samples (values count in-range))
(let ((r (apply 'sample rng distribution-identifier
distribution-parameters)))
(when (< *gslt-lower-limit* r *gslt-upper-limit*)
(multiple-value-bind (bin frac)
(floor (- r *gslt-lower-limit*) *gslt-bin-size*)
(declare (ignore frac))
(incf (aref count bin))
(incf in-range)))))))
(defun total (vector)
(reduce #'+ vector))
(defun chi^2 (pdf-function distribution-identifier
&rest distribution-params)
"Calculate chi-squared for `pdf-function' using
`distribution-params'
Return chi-squred, the normalized cdf, bin counts and total count
modeled after `testpdf' in tests.lisp"
(multiple-value-bind (count/bin num-samples)
(apply #'bin-samples-1 *gslt-bins* *initial-number-of-samples*
distribution-identifier distribution-params)
;; to calculate chi-square, I have to normalize the CDF so that
;; its total equals the number of samples in bins, `num-samples'
(let* ((cdf-norm ; normalized to `num-samples'
(multiple-value-bind (cdf sum)
(distribution-bin-integral-1 pdf-function *gslt-bins*)
(map 'vector (lambda (arg)
(* num-samples
(/ arg sum)))
cdf)))
(chisq (total
(map 'vector (lambda (b c)
(/ (expt (- b c) 2)
c))
count/bin cdf-norm))))
(values chisq num-samples cdf-norm count/bin))))
(defun chi^2-cpdf (chi^2 num-degrees-of-freedom)
"Return the value of the cumulative chi^2 probability function
on `chi^2' (x) with `num-degrees-of-freedom' (a), by calculating
the incomplete gamma function P(a,x)
Here P(a,x) = (incomplete-gamma-function x a)
Values of chi-sq-pdf of 0.95 indicates that the data is drawn from the
same distributions with 95% confidence"
(incomplete-gamma (* 0.5d0 chi^2)
(* 0.5d0 num-degrees-of-freedom)))
;;
(defmacro with-limits ((lower upper) &body body)
`(let* ((*gslt-lower-limit* ,lower)
(*gslt-upper-limit* ,upper)
(*gslt-bin-size* (bin-size *gslt-bins* *gslt-lower-limit*
*gslt-upper-limit*)))
,@body))
;; most tests are done with default limits. For a few this limit is
;; not valid, and I use the `with-limits' macro to specify correct
;; limits
(save-test chisq-tests
(chi^2 (lambda (r) (flat-pdf r *gslt-lower-limit* *gslt-upper-limit*))
:flat :a *gslt-lower-limit* :b *gslt-upper-limit*)
(chi^2 (lambda (r) (tdist-pdf r 1.75d0)) :tdist :nu 1.75d0)
(chi^2 (lambda (r) (tdist-pdf r 12.75d0)) :tdist :nu 12.75d0)
(chi^2 (lambda (r) (cauchy-pdf r 5.0d0)) :levy :c 5.0d0 :alpha 1.0d0) ; levy1
(chi^2 (lambda (r) (cauchy-pdf r 5.0d0)) :levy :c 5.0d0 :alpha 1.01d0) ; levy1a
(chi^2 (lambda (r) (gaussian-pdf r (* (sqrt 2.0d0) 5.0d0)))
:levy :c 5.0d0 :alpha 2.0d0) ; levy2
(chi^2 (lambda (r) (gaussian-pdf r (* (sqrt 2.0d0) 5.0d0)))
:levy :c 5.0d0 :alpha 1.99d0) ; levy2a
(chi^2 (lambda (r) (cauchy-pdf r 5.0d0))
:levy-skew :c 5.0d0 :alpha 1.0d0 :beta 0.0d0) ; levy_skew1
(chi^2 (lambda (r) (cauchy-pdf r 5.0d0))
:levy-skew :c 5.0d0 :alpha 1.01d0 :beta 0.0d0) ; levy_skew1a
(chi^2 (lambda (r) (gaussian-pdf r (* (sqrt 2.0d0) 5.0d0)))
:levy-skew :c 5.0d0 :alpha 2.0d0 :beta 0.0d0) ; levy_skew2
(chi^2 (lambda (r) (gaussian-pdf r (* (sqrt 2.0d0) 5.0d0)))
:levy-skew :c 5.0d0 :alpha 1.99d0 :beta 0.0d0) ; levy_skew2a
(chi^2 (lambda (r) (cauchy-pdf r 5.0d0))
:levy-skew :c 5.0d0 :alpha 1.01d0 :beta 0.001d0) ; levy_skew1b
(chi^2 (lambda (r) (gaussian-pdf r (* (sqrt 2.0d0) 5.0d0)))
:levy-skew :c 5.0d0 :alpha 1.99d0 :beta 0.001d0)
(chi^2 (lambda (r) (laplace-pdf r 2.75d0)) :laplace :a 2.75d0)
(chi^2 (lambda (r) (exponential-power-pdf r 3.7d0 0.3d0)) ; exppow0
:exponential-power :a 3.7d0 :b 0.3d0)
(chi^2 (lambda (r) (exponential-power-pdf r 3.7d0 1.0d0)) ; exppow1
:exponential-power :a 3.7d0 :b 1.0d0)
(chi^2 (lambda (r) (exponential-power-pdf r 3.7d0 1.9d0)) ; exppow1a
:exponential-power :a 3.7d0 :b 1.9d0)
(chi^2 (lambda (r) (exponential-power-pdf r 3.7d0 2.0d0)) ; exppow2
:exponential-power :a 3.7d0 :b 2.0d0)
(chi^2 (lambda (r) (exponential-power-pdf r 3.7d0 3.5d0)) ; exppow2a
:exponential-power :a 3.7d0 :b 3.5d0)
(chi^2 (lambda (r) (exponential-power-pdf r 3.7d0 7.5d0)) ; exppow2b
:exponential-power :a 3.7d0 :b 7.5d0)
(with-limits (0d0 5d0)
(chi^2 (lambda (r) (exponential-pdf r 2.0d0)) :exponential :mu 2.0d0))
(with-limits (0d0 5d0)
(chi^2 (lambda (r) (chisquared-pdf r 13.0d0)) :chi-squared :nu 3.0d0))
(with-limits (0d0 1d0)
(chi^2 (lambda (r) (beta-pdf r 2.0d0 3.0d0)) :beta :a 2.0d0 :b 3.0d0)))
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