;;;; srfi.27.gammas.scm
;;;; Kon Lovett, Dec '17
;;;; Kon Lovett, Jun '17
;;;; Kon Lovett, May '06

; Chicken Generic Arithmetic! (would use fp routines, except for the "real" constraint)

(module srfi.27.gammas

(;export
  *make-random-gammas
  make-random-gammas)

(import scheme
  (chicken base)
  (chicken type)
  (only (type-checks-numbers scheme) check-real)
  (only type-checks-structured check-procedure)
  (srfi 27)
  srfi-27-distributions-support
  (srfi 27 normals))

;;;

(include-relative "srfi-27-common-types")

(: *make-random-gammas (number number random-real-function -> number-function))
(: make-random-gammas (#!rest -> number-function procedure))

;;; Gamma distribution

;NOTE rational vs float "test"
(define-constant *1/3* 1/3)

;; "A Simple Method for Generating Gamma Variables", George Marsaglia & Wai Wan Tsang,
;; ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000, Pages 363 372.

(define (*make-random-gammas alpha theta randoms)
  (if (= 1 alpha)
      ;then special case
      (if (= 1 theta)
          (lambda () (- (log (randoms))))
          (lambda () (* theta (- (log (randoms))))) )
      ;else general case
      (let ((normals (*make-random-normals 0 1 randoms) )
            (uniforms
              (if (< alpha 1)
                  (let ((alpha-inv (*reciprocal alpha)))
                    (lambda () (expt (randoms) alpha-inv) ) )
                  randoms) ) )
        ;
        (let* ((d (- (if (< alpha 1) (+ 1 alpha) alpha) *1/3*) )
               (c (*reciprocal (sqrt (* 9 d))) ) )
          ;
          (lambda ()
            (* theta
               (let loop ()
                 (let* ((x (normals) )
                        (v (+ 1 (* c x)) ) )
                   ;
                   (define (gamma?)
                     (let ((u (uniforms) )
                           (x^2 (* x x) )
                           (v^3 (* v v v) ) )
                       ;
                       (or (< u (- 1 (* 0.0331 x^2 x^2)))
                           (< (log u) (+ (* 0.5 x^2) (* d (- 1 (+ v^3 (log v^3)))))) ) ) )
                   ;
                   (if (and (< 0 v) (gamma?))
                       (* d v)
                       (loop) ) ) ) ) ) ) ) ) )

(define (make-random-gammas #!key (alpha 1) (theta 1) (randoms (current-random-real)))
  (check-positive-real 'make-random-gammas alpha 'alpha)
  (check-positive-real 'make-random-gammas theta 'theta)
  (check-procedure 'make-random-gammas randoms 'randoms)
  (values (*make-random-gammas alpha theta randoms)
          (lambda () (values alpha theta randoms))) )

) ;module srfi.27.gammas