#| 

  Galaxy in Racket Gamble.

  From BLOG example/galaxy.blog
  (run the BLOG model with -s blog.sample.MHSampler to get some result)

var : g0
34994.005105400705: 0.4339999999999998
22929.729460595114: 0.32999999999999996
34638.02450795405: 0.13799999999999996
12117.357091142396: 0.09799999999999995
...
mean: 28721.757314769897

var : g5
14568.675899711576: 0.3929999999999998
11831.433889674547: 0.14299999999999993
14894.598372764063: 0.12399999999999993
...
44085.67436384509: 0.024999999999999988
38504.225550889714: 0.02199999999999999
38911.92059700366: 0.0059999999999999975
mean: 18047.06258982396

var : g9
21805.666821445036: 0.4059999999999998
10340.681733997026: 0.1529999999999999
37162.33526933141: 0.13599999999999993
...
8624.610441950843: 0.03299999999999999
10857.162457725452: 0.02899999999999999
8968.023202696077: 0.0029999999999999988
mean: 24832.632012030208

var : g14
44354.766461949665: 0.2349999999999999
18001.569093746686: 0.21699999999999992
22583.784600120784: 0.2059999999999999
...
39693.14183834324: 0.05099999999999998
36288.14615493976: 0.027999999999999987
34044.41788705041: 0.014999999999999993
mean: 31913.83782650768

var : g15
11217.503478573803: 0.39799999999999985
30249.261083511243: 0.23599999999999993
21813.93255183892: 0.13399999999999995
20138.256062929813: 0.10799999999999994
...
37590.88431343062: 0.09699999999999995
20876.683893695066: 0.02699999999999999
mean: 20911.376860456392

var : g31
47696.77142750669: 0.3669999999999999
14172.368696390813: 0.25699999999999995
43586.689093157496: 0.1809999999999999
20750.694736173496: 0.07099999999999997
...
9873.048460249342: 0.06199999999999998
29313.695207761746: 0.06199999999999998
mean: 32939.0820284139

var : g50
18683.69927192327: 0.5319999999999998
27159.716602000655: 0.3649999999999999
8730.83809879942: 0.07499999999999996
16094.526868001929: 0.016999999999999994
...
43330.33634645631: 0.010999999999999994
mean: 21258.078086370424

var : g70
9388.279285692166: 0.5929999999999999
29710.04477112771: 0.12599999999999995
47144.010147534806: 0.12299999999999996
...
17736.50858648908: 0.027999999999999987
47917.62655062283: 0.02599999999999999
9919.793634516438: 0.014999999999999993
mean: 20154.159962969185

var : g80
29925.028758450873: 0.2739999999999999
48726.28926417515: 0.2579999999999999
40595.73278155002: 0.2099999999999999
45617.45933220525: 0.12799999999999995
...
36202.09901012402: 0.07799999999999996
24232.50804081598: 0.05199999999999998
mean: 39218.833329532594

var : numCluster
14: 1.0
mean: 14.0

var : origCluster0
6: 0.5469999999999999
5: 0.2359999999999999
14: 0.15099999999999994
12: 0.03899999999999998
...
10: 0.02699999999999999
mean: 7.313999999999997


  This is a port of my WebPPL model galaxy.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)
  (; enumerate
   ; rejection-sampler
   ; importance-sampler
   mh-sampler
   
   (define data '(9172 9350 9483 9558 9775 10227 10406 16084 16170 18419 18552 18600 18927 19052 
                       19070 19330 19343 19349 19440 19473 19529 19541 19547 19663 19846 19856 19863 
                       19914 19918 19973 19989 20166 20175 20179 20196 20215 20221 20415 20629 20795 
                       20821 20846 20875 20986 21137 21492 21701 21814 21921 21960 22185 22209 22242 
                       22249 22314 22374 22495 22746 22747 22888 22914 23206 23241 23263 23484 23538 
                       23542 23666 23706 23711 24129 24285 24289 24366 24717 24990 25633 26960 26995 
                       32065 32789 34279))
   
   (define numCluster (add1 (poisson 10.0)))

   (define (clusVelocity c) 
     (uniform 5000.0 50000.0))
   
   
   (define (origCluster g)
     (add1 (random-integer numCluster)))
   
   (define (velocity g) 
     (when (not (= (origCluster g) 0))
       (normal (clusVelocity (origCluster g)) (sqrt 10000))))
   
   (for ([g (length data)])
     (observe-sample (normal-dist (clusVelocity (origCluster g)) (sqrt 10000)) (list-ref data g)))


   (list (clusVelocity (origCluster 0))
         (clusVelocity (origCluster 5))
         (clusVelocity (origCluster 9))
         (clusVelocity (origCluster 14))
         (clusVelocity (origCluster 15))
         (clusVelocity (origCluster 31))
         (clusVelocity (origCluster 50))
         (clusVelocity (origCluster 70))
         (clusVelocity (origCluster 80))
         numCluster
         (origCluster 0)
         )

   )
  )

(show-marginals (model)
                (list  "g0"
                       "g5"
                       "g9"
                       "g14"
                       "g15"
                       "g31"
                       "g50"
                       "g70"
                       "g80"
                       "numCluster"
                       "origCluster0"
                      )
                #:num-samples 1000
                #:truncate-output 3
                ; #:skip-marginals? #t
                ; #:show-stats? #t
                ; #:credible-interval 0.84
                ; #:show-histogram? #t
                ; #:show-percentiles? #t
                )