/* 

  One rigged coin in Picat.

  https://x.com/littmath/status/1827447490683003234
  """
  You have 10000 coins. 9999 of them are fair; one is rigged so that it always 
  lands on heads. You choose a coin at random and flip it 10 times; it’s heads 
  all ten times. The coin is probably
  - fair
  - rigged
  - equaly likely
  - can't tell/see results
  """

  Calculation (for n=10) from
  https://x.com/RRichtsfeld/status/1827476479048880528
  """
  Expected fair all heads: 9999/2^10
  Expected rigged all heads: 1
  Expected total all heads: 9999/2^10+1
  Conditional probability of fair all heads: (9999/2^10)/(9999/2^10+1)≈90.71%
  Conditional probability of rigged all heads: 1/(9999/2^10+1)≈9.29%
  """

  Throwing 10 to 13 heads in a row indicates that fair coin was picked, 
  but when it's 14 heads in a row or more, then we might suspect that it 
  was the biased coin.

  The breakpoint of fair coin -> biased coin can be seen if one compare the 
  probability of selecting a fair coin (1/10000 = 0.0001) with the binomial 
  PDF for throwing n heads in n throws.

  Probability of throwing n=10..13 heads in n flips:

  Picat> X = binomial_dist_pdf(10,1/2,10)
  X = 0.0009765625

  Picat> X = binomial_dist_pdf(11,1/2,11)
  X = 0.00048828125

  Picat> X = binomial_dist_pdf(12,1/2,12)
  X = 0.000244140625

  Picat> X = binomial_dist_pdf(13,1/2,13)
  X = 0.0001220703125

  These is are all greater than the probability of selecting the biased coin (1/10000),
  so we can assume it's probable a fair coin.

  However, the probability of 14 heads in a row is smaller, and we might assume that it's
  actually the biased coin that was selected:

  Picat> X = binomial_dist_pdf(14,1/2,14)
  X = 0.00006103515625

  Pascal Bercker also wrote about this problem when it was published: 
  https://medium.com/@pbercker/given-10000-coins-9999-coins-are-fair-but-1-is-double-headed-fa508fc38d48).
  Using a Netica Bayseian Network model he compared the probability of selecting the 
  biased coin (1/10000) with the probability of getting 10 heads in a rows (1/1024) 
  and (also) concluded the probability of a fair coin is about 90.7%.

  Reading his post, inspired me to do the more detailed analysis of the breakpoint of
  fair/biased coin.

  
  Cf my Gamble model gamble_one_rigged_coin.rkt
  The Gamble model returns the exact probabilities and is much faster than
  this Picat PPL model.

  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.
import util.
% import ordset.

main => go.

/*
  Using observe_abc/2. 

  n = 10
  var : pick coin
  Probabilities:
  fair_coin: 0.9090000000000000
  biased_coin: 0.0910000000000000


  n = 11
  var : pick coin
  Probabilities:
  fair_coin: 0.8420000000000000
  biased_coin: 0.1580000000000000


  n = 12
  var : pick coin
  Probabilities:
  fair_coin: 0.7130000000000000
  biased_coin: 0.2870000000000000


  n = 13
  var : pick coin
  Probabilities:
  fair_coin: 0.5530000000000000
  biased_coin: 0.4470000000000000


  n = 14
  var : pick coin
  Probabilities:
  biased_coin: 0.6150000000000000
  fair_coin: 0.3850000000000000

/*
  Using observe_abc/2. 

  n = 10
  var : pick coin
  Probabilities:
  fair_coin: 0.9090000000000000
  biased_coin: 0.0910000000000000


  n = 11
  var : pick coin
  Probabilities:
  fair_coin: 0.8420000000000000
  biased_coin: 0.1580000000000000


  n = 12
  var : pick coin
  Probabilities:
  fair_coin: 0.7130000000000000
  biased_coin: 0.2870000000000000


  n = 13
  var : pick coin
  Probabilities:
  fair_coin: 0.5530000000000000
  biased_coin: 0.4470000000000000


  n = 14
  var : pick coin
  Probabilities:
  biased_coin: 0.6150000000000000
  fair_coin: 0.3850000000000000

  n = 15
  var : pick coin
  Probabilities:
  biased_coin: 0.7700000000000000
  fair_coin: 0.2300000000000000

  n = 20
  var : pick coin
  Probabilities:
  biased_coin: 1.0000000000000000


*/
go ?=>
  member(N,10..15++[20]),
  println(n=N),
  reset_store,
  run_model(10_000,$model(N),[show_probs,
                              min_accepted_samples=1000 % ,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  fail,
  nl.
go => true.
  
biased_coin(I) = 1. % Only head
fair_coin(I) = categorical([1/2,1/2],[1,0]).

model(N) =>
  Head = 1, Tail = 0,
  % Pick a coin. We have 9999 fair coins and 1 biased coin
  
  PickCoin = cond(flip(9999/10000)==true,fair_coin,biased_coin),
  
  % Throw the coin n times and observe n heads
  Coins = [cond(PickCoin==fair_coin,fair_coin(I),biased_coin(I)) : I in 1..N],
  observe_abc(ones(N,Head),Coins,1,[mean,stdev,sum]),
  
  if observed_ok then
    add("pick coin",PickCoin),
  end.