/* 

  Trick coin in Picat.

  This is a port of (my Gamble port of) the Church model:
  """
  ;; 5.1 The trick coin

  (define observed-data '(h h h h h h h))
  (define num-flips (length observed-data))

  (define samples
    (mh-query
        1000 10

        (define fair-prior 0.999)
        (define fair-coin? (flip fair-prior))

        (define make-coin (lambda (weight) (lambda () (if (flip weight) 'h 't))))
        (define coin (if fair-coin?
                      (make-coin 0.5)
                      (make-coin 0.95)))
                     ;(make-coin (if fair-coin? 0.5 0.95)))
        (define sampled-data (repeat num-flips coin))
        fair-coin?

        (equal? observed-data sampled-data)))

  ;;; (hist samples "Fair coin?")
  """

  Exact probabilities (from my Racket/Gamble model using enumerate):
    (#f : 0.9921875)
    (#t : 0.007812500000000007)


  This is a port of my Racket/Gamble model gamble_trick_coin.rkt

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

*/

import ppl_distributions, ppl_utils.

main => go.

/*

  Var p:
  Probabilities:
  false: 0.99189270343308977 (9910 / 9991)
  true: 0.00810729656691022 (81 / 9991)

*/
go ?=>
  reset_store,
  ObservedData = [h,h,h,h,h,h,h],
  run_model(20_000,$model(ObservedData),[show_probs_rat]),
  nl.
go => true.


/*
  Testing a longer sequence.

  Exact probabilities (from my Gamble model):
  enumerate:
  (#f : 0.99951171875)
  (#t : 0.0004882812500000001)

  This model:
  Var p:
  Probabilities:
  false: 0.99959961563100574 (12483 / 12488)
  true: 0.00040038436899423 (5 / 12488)

*/
go2 ?=>
  reset_store,
  ObservedData = [h,h,h,h,h,h,h,h,h,h,h],
  run_model(100_000,$model(ObservedData),[show_probs_rat,show_histogram]),
  nl.
go2 => true.


make_coin(Weight) = cond(flip(Weight)==true,h,t).
coin(FairCoin) = cond(FairCoin,
                      make_coin(0.5),
                      make_coin(0.95)).
              
model(ObservedData) =>
  NumFlips = ObservedData.len,
  FairPrior = 0.999, 
  FairCoin = flip(FairPrior),
  SampledData = [ coin(FairCoin) : _ in 1..NumFlips],
  P = cond(SampledData==ObservedData,true,false),
  
  observe(FairCoin),
  
  if observed_ok then     
     add("p",P)
   end.