Commit 8abe8ab1 authored by liam's avatar liam
Browse files

New discrete random distributions, make some numbers

double precision.  Remove obsolete vector examples, 
update others.


git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3153 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent 9418a9b8
......@@ -3,7 +3,7 @@
; description: Vectors
; date: Sun Mar 26 2006 - 11:51
; author: Liam M. Healy
; modified: Sat Jun 3 2006 - 21:04
; modified: Sat Nov 11 2006 - 23:20
;********************************************************
;;; $Id: $
......@@ -422,78 +422,48 @@ same length.")
#|
(with-data (vec vector 3)
(with-data (vec vector-double 3)
(setf (gsl-aref vec 0) -3.21d0
(gsl-aref vec 1) 1.0d0
(gsl-aref vec 2) 12.8d0)
(print-vector vec))
(data vec))
(with-data (vec vector 3)
(with-data (vec vector-double 3)
(setf (data vec) #(-3.21d0 1.0d0 12.8d0))
(print-vector vec))
(data vec))
(defparameter vec (make-data 'vector nil 3))
(setf (data vec) #(-3.21d0 1.0d0 12.8d0))
(free vec)
(with-data (vec vector 3)
(with-data (vec vector-double 3)
(setf (data vec) #(-3.21d0 1.0d0 12.8d0))
(data vec))
(with-data (vec vector 3)
(with-data (vec vector-double 3)
(set-all vec 77.8d0)
(print-vector vec)
(set-zero vec)
(print-vector vec))
(data vec))
(with-data (vec vector 3) (set-basis vec 1)
(format t "~&~a ~a ~a"
(gsl-aref vec 0) (gsl-aref vec 1) (gsl-aref vec 2)))
;;; broken
(with-data (vec vector 3)
(with-data (vec vector-double 3)
(setf (gsl-aref vec 0) -3.21d0
(gsl-aref vec 1) 1.0d0
(gsl-aref vec 2) 12.8d0)
(subvector vec 1 2))
(with-data (vec1 vector 3)
(with-data (vec2 vector 3)
(with-data (vec1 vector-double 3)
(with-data (vec2 vector-double 3)
(setf (gsl-aref vec1 0) -3.21d0
(gsl-aref vec1 1) 1.0d0
(gsl-aref vec1 2) 12.8d0)
(vector-copy vec2 vec1)
(print-vector vec2)))
(with-data (vec1 vector 3)
(set-vector vec1 -3.21d0 1.0d0 12.8d0) (print-vector vec1)
(vector-swap-elements vec1 0 1) (print-vector vec1)
(vector-reverse vec1) (print-vector vec1)
(vector*c vec1 2.0d0) (print-vector vec1)
(vector+c vec1 2.0d0) (print-vector vec1))
(with-data (vec1 vector 3)
(with-data (vec2 vector 3)
(set-vector vec1 -3.21d0 1.0d0 12.8d0)
(set-vector vec2 1.4d0 17.0d0 -3.0d0)
(vector+ vec1 vec2)
(print-vector vec1)
(vector* vec1 vec2)
(print-vector vec1)))
(with-data (vec1 vector 3)
(set-vector vec1 -3.21d0 1.0d0 12.8d0)
(print-vector vec1)
(print (vector-min vec1))
(print (vector-max vec1))
(vector-minmax vec1))
(with-data (vec1 vector 3)
(set-vector vec1 -3.21d0 1.0d0 12.8d0)
(print-vector vec1)
(print (vector-min-index vec1))
(print (vector-max-index vec1))
(vector-minmax-index vec1))
(data vec2)))
|#
......@@ -3,7 +3,7 @@
; description: Definition of GSLL system
; date:
; author: Liam Healy
; modified: Sun Oct 29 2006 - 22:17
; modified: Sun Nov 12 2006 - 00:10
;********************************************************
;;; $Id: $
......@@ -118,4 +118,5 @@
(:file "weibull" :depends-on (rng-types))
(:file "gumbel1" :depends-on (rng-types))
(:file "gumbel2" :depends-on (rng-types))
(:file "dirichlet" :depends-on (rng-types))))))
(:file "dirichlet" :depends-on (rng-types))
(:file "discrete" :depends-on (rng-types))))))
......@@ -3,7 +3,7 @@
; description: Dirichlet distribution
; date: Sun Oct 29 2006
; author: Liam M. Healy
; modified: Sun Oct 29 2006 - 22:00
; modified: Sat Nov 11 2006 - 20:57
;********************************************************
;;; $Id: $
......@@ -75,13 +75,13 @@
"0.288000000000d+01"
(with-data (alpha vector-double 4)
(with-data (theta vector-double 4)
(setf (data alpha) #(1.0 2.0 3.0 4.0))
(setf (data theta) #(0.1 0.3 0.4 0.2))
(setf (data alpha) #(1.0d0 2.0d0 3.0d0 4.0d0)
(data theta) #(0.1d0 0.3d0 0.4d0 0.2d0))
(dirichlet-pdf alpha theta))))
(lisp-unit:assert-first-fp-equal
"0.105779029415d+01"
(with-data (alpha vector-double 4)
(with-data (theta vector-double 4)
(setf (data alpha) #(1.0 2.0 3.0 4.0))
(setf (data theta) #(0.1 0.3 0.4 0.2))
(setf (data alpha) #(1.0d0 2.0d0 3.0d0 4.0d0)
(data theta) #(0.1d0 0.3d0 0.4d0 0.2d0))
(dirichlet-log-pdf alpha theta)))))
;********************************************************
; file: discrete.lisp
; description: Discrete random variables
; date: Sat Nov 11 2006 - 21:51
; author: Liam M. Healy
; modified: Sun Nov 12 2006 - 00:09
;********************************************************
;;; $Id: $
(in-package :gsl)
(cffi:defcstruct discrete-t
"Structure for Walker algorithm."
(K :size)
(A :pointer)
(F :pointer))
(defun-gsl discrete-preprocess (probabilities)
"gsl_ran_discrete_preproc"
(((dim0 probabilities) :size) ((gsl-array probabilities) :pointer))
:c-return :pointer
:documentation
"A pointer to a structure that contains the lookup
table for the discrete random number generator. The array probabilities contains
the probabilities of the discrete events; these array elements must all be
positive, but they needn't add up to one (so you can think of them more
generally as ``weights'')---the preprocessor will normalize appropriately.
This return value is used as an argument to #'discrete.")
(defun-gsl discrete-free (table)
"gsl_ran_discrete_free" ((table :pointer))
:c-return :void
:documentation
"De-allocates the lookup table created by #'discrete-preprocess.")
(export 'with-discrete-table)
(defmacro with-discrete-table ((probabilities table) &body body)
`(let ((,table (discrete-preprocess ,probabilities)))
(unwind-protect
(progn ,@body)
(discrete-free ,table))))
(defun-gsl discrete (generator table)
"gsl_ran_discrete"
(((generator generator) :pointer) (table :pointer))
:c-return :size
:documentation
"Generate discrete random numbers after running #'discrete-preprocess;
the argument 'table is the value returned by #'discrete-preprocess.")
(defun-gsl discrete-pdf (k table)
"gsl_ran_discrete_pdf"
((k :size) (table :pointer))
:c-return :double
:documentation
"The probability @math{P[k]} of observing the variable @var{k}.
Since @math{P[k]} is not stored as part of the lookup table, it must be
recomputed; this computation takes @math{O(K)}, so if @var{K} is large
and you care about the original array @math{P[k]} used to create the
lookup table, then you should just keep this original array @math{P[k]}
around.")
;;; Examples and unit test
(lisp-unit:define-test discrete
(lisp-unit:assert-equal
'(1 0 1 1 0 1 1 2 1 2 2)
(with-data (probabilities vector-double 3)
(setf (data probabilities) #(0.25d0 0.5d0 0.25d0))
(with-discrete-table probabilities table
(rng-set *rng-mt19937* 0)
(loop for i from 0 to 10
collect
(discrete *rng-mt19937* table))))
(lisp-unit:assert-first-fp-equal
"0.500000000000d+00"
(with-data (probabilities vector-double 3)
(setf (data probabilities)
#(0.25d0 0.5d0 0.25d0))
(with-discrete-table probabilities
table
(discrete-pdf
1
table))))))
......@@ -3,7 +3,7 @@
; description: Beta distribution
; date: Sun Oct 29 2006
; author: Liam M. Healy
; modified: Sun Oct 29 2006 - 21:07
; modified: Sat Nov 11 2006 - 20:58
;********************************************************
;;; $Id: $
......@@ -70,16 +70,16 @@
(gumbel1 *rng-mt19937* 1.0d0 2.0d0)))))
(lisp-unit:assert-first-fp-equal
"0.296257089650d+00"
(gumbel1-pdf 0.1 1.0 2.0))
(gumbel1-pdf 0.1d0 1.0d0 2.0d0))
(lisp-unit:assert-first-fp-equal
"0.163707359877d+00"
(gumbel1-P 0.1 1.0 2.0))
(gumbel1-P 0.1d0 1.0d0 2.0d0))
(lisp-unit:assert-first-fp-equal
"0.836292640123d+00"
(gumbel1-Q 0.1 1.0 2.0))
(gumbel1-Q 0.1d0 1.0d0 2.0d0))
(lisp-unit:assert-first-fp-equal
"0.100000000000d+00"
(gumbel1-Pinv 0.1637073598773166 1.0 2.0))
(gumbel1-Pinv 0.1637073598773166d0 1.0d0 2.0d0))
(lisp-unit:assert-first-fp-equal
"0.100000000000d+00"
(gumbel1-Qinv 0.8362926401226833 1.0 2.0)))
(gumbel1-Qinv 0.8362926401226833d0 1.0d0 2.0d0)))
......@@ -3,7 +3,7 @@
; description: Beta distribution
; date: Sun Oct 29 2006
; author: Liam M. Healy
; modified: Sun Oct 29 2006 - 21:07
; modified: Sat Nov 11 2006 - 20:58
;********************************************************
;;; $Id: $
......@@ -70,16 +70,16 @@
(gumbel2 *rng-mt19937* 1.0d0 2.0d0)))))
(lisp-unit:assert-first-fp-equal
"0.536256036829d-01"
(gumbel2-pdf 5.0 1.0 2.0))
(gumbel2-pdf 5.0d0 1.0d0 2.0d0))
(lisp-unit:assert-first-fp-equal
"0.818730753078d+00"
(gumbel2-P 10.0 1.0 2.0))
(gumbel2-P 10.0d0 1.0d0 2.0d0))
(lisp-unit:assert-first-fp-equal
"0.181269246922d+00"
(gumbel2-Q 10.0 1.0 2.0))
(gumbel2-Q 10.0d0 1.0d0 2.0d0))
(lisp-unit:assert-first-fp-equal
"1.000000000000d+01"
(gumbel2-Pinv 0.8187307530779818 1.0 2.0))
(gumbel2-Pinv 0.8187307530779818d0 1.0d0 2.0d0))
(lisp-unit:assert-first-fp-equal
"0.100000000000d+02"
(gumbel2-Qinv 0.18126924692201815 1.0 2.0)))
(gumbel2-Qinv 0.18126924692201815d0 1.0d0 2.0d0)))
......@@ -3,7 +3,7 @@
; description: Weibull distribution
; date: Sun Oct 22 2006
; author: Liam M. Healy
; modified: Sun Oct 22 2006 - 22:11
; modified: Sat Nov 11 2006 - 20:59
;********************************************************
;;; $Id: $
......@@ -65,16 +65,16 @@
(weibull *rng-mt19937* 1.0d0 2.0d0)))))
(lisp-unit:assert-first-fp-equal
"0.242631749722d+00"
(weibull-pdf 1.5 1.3 1.0))
(weibull-pdf 1.5d0 1.3d0 1.0d0))
(lisp-unit:assert-first-fp-equal
"0.999288774280d+00"
(weibull-P 3.5 1.3 2.0))
(weibull-P 3.5d0 1.3d0 2.0d0))
(lisp-unit:assert-first-fp-equal
"0.711225720092d-03"
(weibull-Q 3.5 1.3 2.0))
(weibull-Q 3.5d0 1.3d0 2.0d0))
(lisp-unit:assert-first-fp-equal
"0.350000000000d+01"
(weibull-Pinv 0.9992887742799077d0 1.3 2.0))
(weibull-Pinv 0.9992887742799077d0 1.3d0 2.0d0))
(lisp-unit:assert-first-fp-equal
"0.350000000000d+01"
(weibull-Qinv 7.112257200923508e-4 1.3 2.0)))
(weibull-Qinv 7.112257200923508d-4 1.3d0 2.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