;------------------------------------------------------------------------------------------------- ; 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)) ;-------------------------------------------------------------------------------------------------