/* 

  Coin bias in Picat.

  From the R2 model CoinBias.cs 

  Output from the R2 model:
  """
    Mean: 0.421294
    Variance: 0.0162177
    Number of accepted samples = 692
  """

  Cf my Gamble model gamble_coin_bias.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.
import util.

main => go.

/*

  Original sequence: [1,1,0,1,0]

  var : bias
  Probabilities (truncated):
  0.911712951213936: 0.0007800312012480
  0.80827628276366: 0.0007800312012480
  0.803130272025282: 0.0007800312012480
  0.80098415075384: 0.0007800312012480
  .........
  0.092554616366372: 0.0007800312012480
  0.089802945303613: 0.0007800312012480
  0.072189820846068: 0.0007800312012480
  0.051755007423105: 0.0007800312012480
  mean = 0.424589
  variance = 0.0196455
  HPD intervals:
  HPD interval (0.5): 0.31757652819196..0.50724540643249
  HPD interval (0.84): 0.22178253341117..0.62152580556325
  HPD interval (0.9): 0.19478352882288..0.66333135429963
  HPD interval (0.95): 0.16562427568079..0.70465489478389
  HPD interval (0.99): 0.08980294530361..0.75272087129763
  HPD interval (0.999): 0.05175500742310..0.80827628276366
  HPD interval (0.9999): 0.05175500742310..0.91171295121394
  HPD interval (0.99999): 0.05175500742310..0.91171295121394



  Longer sequence (bern_n(0.12,20) ->  [0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0])

  var : bias
  Probabilities (truncated):
  0.3805668223137: 0.0024875621890547
  0.29698129706488: 0.0024875621890547
  0.296227731886613: 0.0024875621890547
  0.279055354069077: 0.0024875621890547
  .........
  0.013035569418133: 0.0024875621890547
  0.012786837453995: 0.0024875621890547
  0.011346665651321: 0.0024875621890547
  0.007748917765183: 0.0024875621890547
  mean = 0.112944
  variance = 0.0035595
  HPD intervals:
  HPD interval (0.5): 0.04872648794913..0.12220696207965
  HPD interval (0.84): 0.02305950387306..0.17467245843439
  HPD interval (0.9): 0.01557719650712..0.20235048741404
  HPD interval (0.95): 0.01557719650712..0.23057978553707
  HPD interval (0.99): 0.00774891776518..0.27241727394527
  HPD interval (0.999): 0.00774891776518..0.38056682231370
  HPD interval (0.9999): 0.00774891776518..0.38056682231370
  HPD interval (0.99999): 0.00774891776518..0.38056682231370

*/
go ?=>
  member(X,[[1,1,0,1,0],[0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]]),
  println(x=X),
  reset_store,
  run_model(100_000,$model(X),[show_probs_trunc,mean,variance,show_hpd_intervals]),
  nl,
  fail,
  nl.
go => true.
  
  
model(X) =>
  % X = [1,1,0,1,0],
  
  % Longer sequence
  % X=bern_n(0.12,20)
  % X = [0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

  
  % beta_dist(2,5) has mean about 0.2855
  % Bias = beta_dist(2,5),
  Bias = beta_dist(2,5),
  foreach(V in X)
    observe(bern(Bias)==V)  
  end,
  
  if observed_ok then
    add("bias",Bias),
  end.