source: project/release/4/srfi-27/trunk/srfi-27-gammas.scm @ 35456

Last change on this file since 35456 was 35456, checked in by Kon Lovett, 17 months ago

begin typing, dep /

File size: 2.4 KB
Line 
1;;;; srfi-27-gammas.scm
2;;;; Kon Lovett, Dec '17
3;;;; Kon Lovett, Jun '17
4;;;; Kon Lovett, May '06
5
6; Chicken Generic Arithmetic! (could use fp routines)
7
8(module srfi-27-gammas
9
10(;export
11  *make-random-gammas
12  make-random-gammas)
13
14(import scheme chicken)
15
16(use
17  (only type-errors error-argument-type)
18  (only type-checks
19    define-check+error-type
20    check-procedure
21    check-cardinal-integer
22    check-real
23    check-open-interval
24    check-closed-interval)
25  srfi-27
26  srfi-27-distributions-support
27  srfi-27-normals)
28
29;;; Gamma distribution
30
31;; "A Simple Method for Generating Gamma Variables", George Marsaglia & Wai Wan Tsang,
32;; ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000, Pages 363 372.
33
34(define (*make-random-gammas alpha theta randoms)
35  (if (= 1.0 alpha)
36    ;then special case
37    (if (= 1.0 theta)
38      (lambda () (- (log (randoms))))
39      (lambda () (* theta (- (log (randoms))))) )
40    ;else general case
41    (let (
42      (normals (*make-random-normals 0.0 1.0 randoms) )
43      (uniforms
44        (if (< alpha 1.0)
45          (let ((alpha-inv (*reciprocal alpha)))
46            (lambda () (expt (randoms) alpha-inv) ) )
47          randoms) ) )
48      ;
49      (let* (
50        (d (- (if (< alpha 1.0) (+ 1.0 alpha) alpha) *one-third*) )
51        (c (*reciprocal (sqrt (* 9.0 d))) ) )
52        ;
53        (lambda ()
54          (*
55            theta
56            (let loop ()
57              (let* (
58                (x (normals) )
59                (v (+ 1.0 (* c x)) ) )
60                ;
61                (define (gamma?)
62                  (let (
63                    (u (uniforms) )
64                    (x^2 (* x x) )
65                    (v^3 (* v v v) ) )
66                    ;
67                    (or
68                      (< u (- 1.0 (* 0.0331 x^2 x^2)))
69                      (< (log u) (+ (* 0.5 x^2) (* d (- 1.0 (+ v^3 (log v^3)))))) ) ) )
70                ;
71                (if (and (< 0.0 v) (gamma?))
72                   (* d v)
73                   (loop) ) ) ) ) ) ) ) ) )
74
75(define (make-random-gammas #!key (alpha 1.0) (theta 1.0) (randoms (current-random-real)))
76  (check-positive-real 'make-random-gammas alpha 'alpha)
77  (check-positive-real 'make-random-gammas theta 'theta)
78  (check-procedure 'make-random-gammas randoms 'randoms)
79  (values
80    (*make-random-gammas alpha theta randoms)
81    (lambda () (values alpha theta randoms))) )
82
83) ;module srfi-27-gammas
Note: See TracBrowser for help on using the repository browser.