;-------------------------------------------------------------------------------------------------
; PCM, my_gamma --- Gamma Distribution Gamma(a, b) = Gamma(shape, 1/rate) ---
;                   mean = shape/rate = a/(1/b) = a*b
;                   std  = sqrt(shape)/rate = sqrt(a)/(1/b) = sqrt(a)*b
;                   Fig.9.8 is in Kruschke, Doing Bayesian Data Analysis, 2011, p.209
; PCM 2014/08/08
;-------------------------------------------------------------------------------------------------
(define n 60000)
(define (my_gamma a b string)
  (let* ((comment1 (display "Gamma(a , b) = Gamma(shape, 1/rate) = " (pair a b)))
         (expected-mean (* a b))
         (comment2 (display "expected mean = " expected-mean))
         (std (* (sqrt a) b))
         (comment3 (display "expected std = " std))
         (comment4 (display "n of samples = " n))
         (vals (repeat n (lambda () (gamma a b))))
         (sample-based-mean (mean vals))
         (comment5 (display "sample-based mean = " sample-based-mean))
         (deviation (abs (- expected-mean sample-based-mean)))
         (comment6 (display "|deviation of means| = " deviation))
         (comment7 (density vals string))
         )
    (display "################################################################################")))

(define start_time (get-time))
(my_gamma 0.5  (/ 1.0 0.0044) "G(shape,rate) = G(0.5, 0.0044) - Fig.9.8a, Krusche,2011")
(my_gamma 0.51 (/ 1.0 0.01  ) "G(shape,rate) = G(0.51,0.01)   - Fig.9.8b, Krusche,2011") 
(my_gamma 1.0  (/ 1.0 0.1   ) "G(shape,rate) = G(1.0, 0.1)    - Fig.9.8c, Krusche,2011") 
(my_gamma 1.2  (/ 1.0 0.025 ) "G(shape,rate) = G(1.2, 0.025)  - Fig.9.8d, Krusche,2011") 
(my_gamma 1.0  (/ 1.0 1.0   ) "G(shape,rate) = G(1.0, 1.0)    - Fig.9.8e, Krusche,2011") 
(my_gamma 6.2  (/ 1.0 0.12  ) "G(shape,rate) = G(6.2, 0.12)   - Fig.9.8f, Krusche,2011") 
(display "computation time in sec =" (/ (- (get-time) start_time) 1000))
;-------------------------------------------------------------------------------------------------