/* 

  Cheating in Picat.

  From 
  https://discourse.julialang.org/t/fitting-a-observed-value-to-a-binomial-distribution-turing/66619
  """
  Fitting a observed value to a Binomial Distribution Turing

  I am new to Turing and trying to learn it by trying to replicate the Chapters from the book 
  https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers 2 .

  Here is my problem (from Chapter 2) :

  1.) Have a probability p ~ Uniform(0,1) which I have to infer the “number_of_cheaters”
  2.) p_skewed = p*(0.5) + (0.25) (deterministic value)
  3.) with model:
  yes_responses = pm.Binomial(“number_of_cheaters”, 100, p_skewed, observed=35)

  How do I write this in Turing ?
  """

  Cf my Gamble model gamble_cheating2.rkt
  Also, cf ppl_cheating.pi

  This is a port of my WebPPL model cheating2.wppl (via the Gamble model)
  which is a port of an adaped Turing.jl model that I wrote as the answer.

  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.

/*

  var : p
  Probabilities (truncated):
  0.463000492874068: 0.0047619047619048
  0.45652546242649: 0.0047619047619048
  0.424643612198831: 0.0047619047619048
  0.419178425063928: 0.0047619047619048
  .........
  0.052477955376952: 0.0047619047619048
  0.038879434596225: 0.0047619047619048
  0.031304243966613: 0.0047619047619048
  0.023448659583669: 0.0047619047619048
  mean = 0.208359

  var : p skewed
  Probabilities (truncated):
  0.481500246437034: 0.0047619047619048
  0.478262731213245: 0.0047619047619048
  0.462321806099416: 0.0047619047619048
  0.459589212531964: 0.0047619047619048
  .........
  0.276238977688476: 0.0047619047619048
  0.269439717298113: 0.0047619047619048
  0.265652121983307: 0.0047619047619048
  0.261724329791835: 0.0047619047619048
  mean = 0.35418


*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean]),
  nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  N = 100,
  YesResponses = 35,

  P = uniform(0,1),
  PSkewed = (P * 0.5) + 0.25,

  observe(binomial_dist(N,PSkewed)==YesResponses),
  
  
  if observed_ok then
    add("p",P),
    add("p skewed",PSkewed),    
  end.