#| 

  Gaussian mixture model in Racket Gamble.

var : a
0.15779851745397125: 0.0022382282728588493
0.8222299073878259: 0.0022356681420079806
0.5821844600360766: 0.002213696264596593
0.4399328626473517: 0.0022099641622077976
...
-0.6318480927553968: 1.5843953294680266e-5
-0.39677567093850574: 1.580480413806597e-5
-0.29334073071714306: 1.5463267024338e-5
-0.5162286882697529: 1.4617786641442041e-5
mean: 0.15223049496434465

var : b
0.7784775104195165: 0.0022382282728588493
0.31184643224441866: 0.0022356681420079806
0.8292299279663305: 0.002213696264596593
0.9491608069337552: 0.0022099641622077976
...
-0.9849728382831324: 1.5843953294680266e-5
-0.9853692657679337: 1.580480413806597e-5
-0.9888655796377084: 1.5463267024338e-5
-0.9978285656656004: 1.4617786641442041e-5
mean: 0.3313984343036442

var : p
0.3907427976129233: 0.0022382282728588493
0.0774098697523533: 0.0022356681420079806
0.56749734976977: 0.002213696264596593
0.15181248402782035: 0.0022099641622077976
...
0.010953324593766583: 1.5843953294680266e-5
0.07947177684685489: 1.580480413806597e-5
0.14870647317126526: 1.5463267024338e-5
0.08454344524640289: 1.4617786641442041e-5
mean: 0.3159718935367733

var : a-b
-0.6206789929655453: 0.0022382282728588493
0.5103834751434072: 0.0022356681420079806
-0.24704546793025395: 0.002213696264596593
-0.5092279442864034: 0.0022099641622077976
...
0.35312474552773565: 1.5843953294680266e-5
0.588593594829428: 1.580480413806597e-5
0.6955248489205653: 1.5463267024338e-5
0.48159987739584753: 1.4617786641442041e-5
mean: -0.1791679393392997

var : a>b
#f: 0.5793938740122834
#t: 0.42060612598771685
mean: 0.42060612598771685

  This is a port of my WebPPL model gaussian_mixture_model.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 (gaussian-mixture-problem)
  (; enumerate
   ; rejection-sampler
   importance-sampler
   ; mh-sampler

   (define xs '(0.2 1.0 0.5 0.6))
    
   (define p (beta 0.5 1))
   (define (z i) (bernoulli-dist p))
    
   (define a (uniform -1 1))
   (define b (uniform -1 1))   
   
   (define (x i) 
     (if (= (sample (z i)) 1)
         (normal-dist a 1.0)
         (normal-dist b 1.0)))
    
   (observe-sample (x 0) 0.2)
   (observe-sample (x 1) 1.0)   
   (observe-sample (x 3) 0.5)
   (observe-sample (x 4) 0.6)
   
   (list a b p (- a b) (> a b))
   )
  )

(show-marginals (gaussian-mixture-problem)
                (list "a"
                      "b"
                      "p"
                      "a-b"
                      "a>b"
                      )
                #:num-samples 1000
                #:truncate-output 4
                )