#| Binomial Trial Count in Racket Gamble. From infer.net test/Tests/BlogTests.cs """ Example of inferring the size of a population. Reference: "A. Raftery, "Inference for the binomial N parameter: A hierarchical Bayes approach", Biometrika 1988 http://pluto.huji.ac.il/~galelidan/52558/Material/Raftery.pdf .... ExpectationPropagation impala theta = Beta(80.24,91.17)(mean=0.4681) N mode = 44, median = 45 waterbuck theta = Beta(238,233.8)(mean=0.5044) N mode = 124, median = 125 impala N mode = 37, median = 65 waterbuck N mode = 122, median = 208 VariationalMessagePassing impala theta = Beta(106,370.9)(mean=0.2223) N mode = 94, median = 94 waterbuck theta = Beta(316,393)(mean=0.4457) N mode = 141, median = 141 impala N mode = 37, median = 65 waterbuck N mode = 122, median = 208 """ * Test: impala var : N 62: 0.05500816241399386 53: 0.05440910839260953 33: 0.04826823738023672 42: 0.047442772091941476 57: 0.04571031765479186 ... 385: 1.5740013030021933e-298 268: 3.6923792390778295e-304 363: 4.90937576822e-312 326: 1.866409666e-315 320: 1.27010167e-316 mean: 84.86887557077978 Credible interval (0.84): 27..104 var : theta 0.3450645496275337: 0.02989459985923875 0.4992095322255353: 0.02957658001386942 0.36432762206475977: 0.029510379242857303 0.3861226192419058: 0.029281248591698023 0.3581017286781968: 0.028970083373482506 ... 0.549381947753365: 7.388558191e-313 0.780905513977937: 6.7189243348e-313 0.49314641442548296: 1.866409666e-315 0.5007740892785166: 1.27010167e-316 0.43660500087850623: 1.169394e-316 mean: 0.4177894719709287 Credible interval (0.84): 0.16962914225347447..0.621518794702926 var : beta_1 29.75567868826473: 0.02989459985923875 68.8950734665327: 0.02957658001386942 23.684869148035723: 0.029510379242857303 54.172711211425245: 0.029281248591698023 24.000091999633973: 0.028970083373482506 ... 52.13708460012303: 7.388558191e-313 96.24434227396809: 6.7189243348e-313 8.96239894427801: 1.866409666e-315 80.95175000798517: 1.27010167e-316 92.02566606154166: 1.169394e-316 mean: 46.521818094282565 Credible interval (0.84): 1.9745322656824051..75.36690477603493 var : beta_2 51.248827701098335: 0.02989459985923875 77.61548225488987: 0.02957658001386942 52.38787876152872: 0.029510379242857303 95.07917019875428: 0.029281248591698023 36.09802502193238: 0.028970083373482506 ... 31.95423102741597: 7.388558191e-313 29.84262273154817: 6.7189243348e-313 14.743251734389077: 1.866409666e-315 82.23714285051119: 1.27010167e-316 95.67703946019157: 1.169394e-316 mean: 58.001376342452836 Credible interval (0.84): 25.63973921359651..92.6405239656868 * Test: waterbuck var : N 313: 0.0658060500995476 118: 0.047164229370357386 251: 0.043047513696568965 209: 0.0425769531713656 574: 0.04204287279131081 ... 368: 1.5735326637442952e-300 487: 1.3373606693786996e-300 768: 9.052385967269911e-304 449: 1.289029268815856e-304 590: 1.205444231251913e-309 mean: 211.63834624572544 Credible interval (0.84): 80..407 var : theta 0.20321987938598077: 0.04341075307794026 0.2539473152262796: 0.043047513696568965 0.2979177756885858: 0.0425769011286651 0.11106753761834746: 0.04204287279131081 0.3707602848809138: 0.041899013409286694 ... 0.33533588357591276: 9.052385967269911e-304 0.5154217186731817: 1.289029268815856e-304 0.8249148321804574: 4.567248240797582e-306 0.46225246083677096: 1.9321799176946618e-307 0.42046072514365096: 1.205444231251913e-309 mean: 0.3760298011302877 Credible interval (0.84): 0.11155196723654083..0.5534364675050459 var : beta_1 24.57322278449553: 0.04341075307794026 27.04533735567941: 0.043047513696568965 25.567059162107373: 0.0425769011286651 8.486855972177825: 0.04204287279131081 71.09657088955555: 0.041899013409286694 ... 15.38376331776911: 9.052385967269911e-304 98.601441407809: 1.289029268815856e-304 99.74613536123563: 4.567248240797582e-306 60.19387679326497: 1.9321799176946618e-307 16.71942892364255: 1.205444231251913e-309 mean: 41.82330689783282 Credible interval (0.84): 4.959553625664477..77.35937930400272 var : beta_2 53.59955035166501: 0.04341075307794026 96.24383479203509: 0.043047513696568965 42.21999343018016: 0.0425769011286651 41.932582253910866: 0.04204287279131081 97.36974188897442: 0.041899013409286694 ... 31.4627689671833: 9.052385967269911e-304 94.82733260350889: 1.289029268815856e-304 27.847440462528645: 4.567248240797582e-306 87.68285213444226: 1.9321799176946618e-307 45.973264935132846: 1.205444231251913e-309 mean: 64.63182967307834 Credible interval (0.84): 30.97885358987878..99.14916313656745 This is a port of my WebPPL model binomial_trial_count.wppl 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") (define (model test) (; enumerate ; rejection-sampler importance-sampler ; mh-sampler (define maxN 1000) ; random Boolean evidence ~ BooleanDistrib(0.5); (define beta_1 (uniform 0.1 100)) (define beta_2 (uniform 0.1 100)) (define theta (beta beta_1 beta_2)) (define priors (for/list ([i maxN]) (cons i (if (= i 0) 0 (/ 1 i))))) (define N (add1 (discrete priors))) ;; (define N (add1 (random-integer maxN))) ;; Faster but strange results (define (x i) (binomial-dist N theta)) (define impala '(15 20 21 23 26)) (define waterbuck '(53 57 66 67 72)) ;; Impala (if (eq? test "impala") (for ([i (length impala)]) (observe-sample (x i) (list-ref impala i)) ) (for ([i (length waterbuck)]) (observe-sample (x i) (list-ref waterbuck i)) ) ) (list N theta beta_1 beta_2 ) ) ) (for ([test '("impala" "waterbuck")]) (show " * Test" test) (show-marginals (model test) (list "N" "theta" "beta_1" "beta_2" ) #:num-samples 1000 #:truncate-output 5 ; #:skip-marginals? #t ; #:show-stats? #t #:credible-interval 0.84 ; #:show-histogram? #t ; #:show-percentiles? #t ) )