#| 

  Poisson dist recover parameters in Racket/Gamble 

  Generating random (poisson 20) and trying to recover data

data length:: 250
(min: 8 mean: 19.456 max: 33)
variable : lambda
19.45650601844706: 0.0038014941387805
19.45540257970331: 0.0038014941387805
19.455863742677668: 0.0038014941387805
19.45959412114821: 0.0038011726457669454
19.45208930017503: 0.0038010118992601676
...
20.77217382306073: 1.6074650677747472e-7
18.24097154436201: 1.6074650677747472e-7
20.750043003328646: 1.6074650677747472e-7
20.726426775934648: 1.6074650677747472e-7
15.368817536249583: 0.0
mean: 19.45011621405835
Min: 18.497331726843033 Mean: 19.44668087997517 Max: 20.343966790456953 Variance: 0.05438304837292746 Stddev: 0.23320173321167117
HPD interval (0.84): 19.12837193411687..19.75441227899537
HPD interval (0.9): 19.08445108298639..19.841138387795684
HPD interval (0.95): 19.050308051523363..19.924805789843433
HPD interval (0.99): 18.915593911109504..20.036776727318237
Percentiles:
(0.01 19.00263244376388)
(0.025 19.02459729906923)
(0.1 19.162279189510752)
(0.05 19.10032442696338)
(0.25 19.28019131946778)
(0.5 19.44474915295976)
(0.75 19.59807770672182)
(0.84 19.67785876845436)
(0.9 19.75114276618536)
(0.95 19.863535253070793)
(0.975 19.905662795296408)
(0.99 19.97305019037114)
(0.999 20.186241540973203)
Histogram:
18.497:   1 # (0.001 / 0    )
18.59 :   0  (0     / 0.001)
18.682:   2 # (0.002 / 0.001)
18.774:   0  (0     / 0.003)
18.867:   3 ## (0.003 / 0.003)
18.959:   3 ## (0.003 / 0.006)
19.051:  24 ########### (0.024 / 0.009)
19.144:  38 ################## (0.038 / 0.033)
19.236: 125 ########################################################## (0.125 / 0.071)
19.328: 129 ########################################################### (0.129 / 0.196)
19.421: 134 ############################################################## (0.134 / 0.325)
19.513: 175 ################################################################################ (0.175 / 0.459)
19.605: 121 ######################################################## (0.121 / 0.634)
19.698: 102 ############################################### (0.102 / 0.755)
19.79 :  61 ############################ (0.061 / 0.857)
19.882:  37 ################# (0.037 / 0.918)
19.975:  35 ################ (0.035 / 0.955)
20.067:   7 #### (0.007 / 0.99 )
20.159:   1 # (0.001 / 0.997)
20.252:   1 # (0.001 / 0.998)

variable : post
20: 0.0970265914908837
21: 0.09408718086795081
18: 0.08549929874336426
19: 0.08007940877434813
17: 0.07248494006164635
...
9: 0.002953877808542876
6: 0.0024486515377412726
40: 0.0007563123143880186
34: 0.00020125462648539836
33: 6.445934921776736e-5
mean: 19.83709080972036
Min: 3 Mean: 19.384 Max: 32 Variance: 19.546544 Stddev: 4.421147362393613
HPD interval (0.84): 14..25
HPD interval (0.9): 11..25
HPD interval (0.95): 10..27
HPD interval (0.99): 9..30
Percentiles:
(0.01 10)
(0.025 11)
(0.1 14)
(0.05 13)
(0.25 16)
(0.5 19)
(0.75 23)
(0.84 24)
(0.9 25)
(0.95 27)
(0.975 28)
(0.99 30)
(0.999 32)
Histogram:
 3:  8 ####### (0.008 / 0    )
 9:  1 # (0.001 / 0.008)
10:  8 ####### (0.008 / 0.009)
11: 12 ########### (0.012 / 0.017)
12: 12 ########### (0.012 / 0.029)
13: 15 ############# (0.015 / 0.041)
14: 34 ############################# (0.034 / 0.056)
15: 52 ############################################# (0.052 / 0.09 )
16: 65 ######################################################## (0.065 / 0.142)
17: 64 ####################################################### (0.064 / 0.207)
18: 94 ################################################################################ (0.094 / 0.271)
19: 84 ######################################################################## (0.084 / 0.365)
20: 79 #################################################################### (0.079 / 0.449)
21: 77 ################################################################## (0.077 / 0.528)
22: 87 ########################################################################### (0.087 / 0.605)
23: 65 ######################################################## (0.065 / 0.692)
24: 58 ################################################## (0.058 / 0.757)
25: 76 ################################################################# (0.076 / 0.815)
26: 44 ###################################### (0.044 / 0.891)
27: 21 ################## (0.021 / 0.935)
28: 14 ############ (0.014 / 0.956)
29: 16 ############## (0.016 / 0.97 )
30:  6 ###### (0.006 / 0.986)
32:  8 ####### (0.008 / 0.992)

variable : p
#f: 0.5684679652305308
#t: 0.43153203476946955
mean: 0.43153203476946955
Min: 0 Mean: 0.403 Max: 1 Variance: 0.240591 Stddev: 0.49050076452539804
Percentiles:
(0.01 #f)
(0.025 #f)
(0.1 #f)
(0.05 #f)
(0.25 #f)
(0.5 #f)
(0.75 #t)
(0.84 #t)
(0.9 #t)
(0.95 #t)
(0.975 #t)
(0.99 #t)
(0.999 #t)
Histogram:
#f: 597 ################################################################################ (0.597 / 0    )
#t: 403 ####################################################### (0.403 / 0.597)


  This program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Racket page: http://www.hakank.org/racket/

|#

#lang gamble

; (require gamble/viz)
(require racket)
(require "gamble_utils.rkt")
; (require "gamble_distributions.rkt")

; Note: For more than about 250, importance-sampler returns a lot of +nan.0
(define num-samples 250)
(define data (for/list ((i num-samples)) (poisson 20)))

(show "data length:" (length data))
(show2 "min:" (apply min data) "mean:" (* 1.0 (avg data)) "max:" (apply max data))

(define (model)
  (;enumerate
   ; rejection-sampler
   importance-sampler
   ; mh-sampler

   ; (define lambda_ (add1 (random-integer 100)))
   (define lambda_ (normal (avg data) 1))
    
   (for ([i num-samples])
     (observe-sample (poisson-dist lambda_) (list-ref data i)))

   (define post (poisson lambda_))
   (define p (> post 20))
    
   (list lambda_
         post
         p
    )
    
   )
)

(show-marginals (model)
                (list  "lambda"
                       "post"
                       "p"
                     )
                    #:num-samples 1000
                    #:truncate-output 5
                    ; #:skip-marginals? #t
                    #:show-stats? #t
                    ; #:credible-interval 0.84
                    #:hpd-interval (list 0.84 0.9 0.95 0.99)
                    #:show-histogram? #t
                    #:show-percentiles? #t
                    ; #:burn 0
                    ; #:thin 0
                    )