#| 

  Loaded coin in Racket Gamble.

  From Osvaldo Martin: "Bayesian Analysis with Python", 2nd edition,
  page 39, exercise 6.
  """
  Let's suppose that we have two coins; when we toss the first coint, half
  of the time it lands tails and half of the time on heads. The other coin
  is a loaded coin that always lands on heads. If we take one if the coins
  at random and get a head, what is the probability the the coin is the
  unfair one?
  """    

  Model 1
  var : coin = fair
  #f: 2/3 (0.6666666666666666)
  #t: 1/3 (0.3333333333333333)
  mean: 1/3 (0.3333333333333333)

  var : coin = loaded
  #t: 2/3 (0.6666666666666666)
  #f: 1/3 (0.3333333333333333)
  mean: 2/3 (0.6666666666666666)


  Model 2
  var : coin = fair
  #f: 2/3 (0.6666666666666666)
  #t: 1/3 (0.3333333333333333)
  mean: 1/3 (0.3333333333333333)

  var : coin = loaded
  #t: 2/3 (0.6666666666666666)
  #f: 1/3 (0.3333333333333333)
  mean: 2/3 (0.6666666666666666)


  This is a port of my PSI model loaded_coin.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")

; This model uses set!.
; See below (model2) for a set!-free model
(define (model)
  
  (enumerate ; cannot enumerate irichlet-dist
   ; rejection-sampler
   ; importance-sampler
   ; mh-sampler

   ;; What coin?
   (define fair 0)
   (define loaded 1)

   ; The outcome
   (define tail 0)
   (define head 1)

   ; coin := [fair,loaded][categorical([1/2,1/2])];
   (define coin (categorical-vw2 (vector 1/2 1/2) (vector fair loaded)))
   
   (define outcome 0)
   (if (= coin fair)
       (set! outcome (categorical-vw2 (vector 1/2 1/2) (vector tail head))) ; [tail,head][categorical([1/2,1/2])];
       (set! outcome head))

   (observe/fail (= outcome head))
   
   (list (= coin fair)
         (= coin loaded)
         )

   
   )
  )

(displayln "Model 1")
(show-marginals (model)
                (list "coin = fair"
                      "coin = loaded"
                      )
                #:num-samples 1000
                ; #:truncate-output 1
                ; #:skip-marginals? #t
                ; #:credible-interval 0.94
                )

; Without set!
(define (model2)
  
  (enumerate ; cannot enumerate irichlet-dist
   ; rejection-sampler
   ; importance-sampler
   ; mh-sampler

   ;; What coin?
   (define fair 0)
   (define loaded 1)

   ; The outcome
   (define tail 0)
   (define head 1)

   ; coin := [fair,loaded][categorical([1/2,1/2])];
   (define coin (categorical-vw2 (vector 1/2 1/2) (vector fair loaded)))
   
   (define outcome (if (= coin fair)
                       (categorical-vw2 (vector 1/2 1/2) (vector tail head))
                       head))
     
   (observe/fail (= outcome head))
   
   (list (= coin fair)
         (= coin loaded)
         )
   
   )
  )

(displayln "\nModel 2")
(show-marginals (model2)
                (list "coin = fair"
                      "coin = loaded"
                      )
                #:num-samples 1000
                ; #:truncate-output 1
                ; #:skip-marginals? #t
                ; #:credible-interval 0.94
                )