/* 

  Handedness in Picat.

  From infer.net test/Tests/BlogTests.cs
  """
  probRightHandedActual: Beta(7.72,3.08)[mean=0.7148]
  """

  This is a port of my Gamble model gamble_handesness.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.

/*

data = [false,true,true,true,true,true,true,true,false,false]
var : betaA
Probabilities (truncated):
9.733654460745703: 0.0138888888888889
9.558978215585919: 0.0138888888888889
9.298318735928424: 0.0138888888888889
9.229999463646672: 0.0138888888888889
.........
0.856503765497591: 0.0138888888888889
0.72743710164327: 0.0138888888888889
0.718599609433953: 0.0138888888888889
0.538810608228115: 0.0138888888888889
mean = 5.76192

var : betaB
Probabilities (truncated):
9.166293255596559: 0.0138888888888889
8.99673291435313: 0.0138888888888889
8.346112579268455: 0.0138888888888889
7.547421938528968: 0.0138888888888889
.........
0.732565033590684: 0.0138888888888889
0.561238266789: 0.0138888888888889
0.482516978160719: 0.0138888888888889
0.331073366259725: 0.0138888888888889
mean = 3.89065

var : prob
Probabilities (truncated):
0.946064656052852: 0.0138888888888889
0.924345165588011: 0.0138888888888889
0.903690913363414: 0.0138888888888889
0.878286383572196: 0.0138888888888889
.........
0.417321538626373: 0.0138888888888889
0.411852922857473: 0.0138888888888889
0.404101423093683: 0.0138888888888889
0.396885375774647: 0.0138888888888889
mean = 0.639645

data = [false,true,true,true,true,true,true,true,true,true]
var : betaA
Probabilities (truncated):
9.997325581497199: 0.0014306151645207
9.981977380803777: 0.0014306151645207
9.977994700883512: 0.0014306151645207
9.977803705249821: 0.0014306151645207
.........
0.540053001856456: 0.0014306151645207
0.311752115521465: 0.0014306151645207
0.306443558217233: 0.0014306151645207
0.067098215253604: 0.0014306151645207
mean = 6.37659

var : betaB
Probabilities (truncated):
9.733448047066782: 0.0014306151645207
9.642609571871631: 0.0014306151645207
9.636759930121135: 0.0014306151645207
9.49398410017322: 0.0014306151645207
.........
0.096171544909557: 0.0014306151645207
0.085036936255655: 0.0014306151645207
0.083257504777637: 0.0014306151645207
0.045142229667465: 0.0014306151645207
mean = 2.76029

var : prob
Probabilities (truncated):
0.996708490349527: 0.0014306151645207
0.994738072508563: 0.0014306151645207
0.994584869604441: 0.0014306151645207
0.993363980143458: 0.0014306151645207
.........
0.397959757730835: 0.0014306151645207
0.339175260029788: 0.0014306151645207
0.32298410979033: 0.0014306151645207
0.321111906743192: 0.0014306151645207
mean = 0.813889

data = [false,true,true,true,true,true,true,true,true,true,false,true,true,true,true,true,true,true,true,true]
var : betaA
Probabilities (truncated):
19.8785956482769: 0.0526315789473684
17.405564634783921: 0.0526315789473684
17.276545184327542: 0.0526315789473684
17.230718618831933: 0.0526315789473684
.........
8.377147944819717: 0.0526315789473684
8.213454544643618: 0.0526315789473684
6.633435556028707: 0.0526315789473684
2.867108025991874: 0.0526315789473684
mean = 12.8203

var : betaB
Probabilities (truncated):
9.399792174529187: 0.0526315789473684
8.293329881640771: 0.0526315789473684
8.015176173306619: 0.0526315789473684
7.155302552112984: 0.0526315789473684
.........
0.93469064726247: 0.0526315789473684
0.868122820215357: 0.0526315789473684
0.852891859064294: 0.0526315789473684
0.825163554784452: 0.0526315789473684
mean = 3.76133

var : prob
Probabilities (truncated):
0.924953216152386: 0.0526315789473684
0.919385545282276: 0.0526315789473684
0.90467062934394: 0.0526315789473684
0.904581019599599: 0.0526315789473684
.........
0.813442636622623: 0.0526315789473684
0.763098799249227: 0.0526315789473684
0.73955430972024: 0.0526315789473684
0.672312827186952: 0.0526315789473684
mean = 0.854876

*/
go ?=>
  StudentData = [false,true,true,true,true,true,true,true,false,false],
  LecturerData = [false,true,true,true,true,true,true,true,true,true],
  % A longer sequences
  TestData = LecturerData ++ LecturerData,  
  member(Data,[StudentData,LecturerData,TestData]),
  println(data=Data),
  reset_store,
  run_model(100_000,$model(Data),[show_probs_trunc,mean]),
  fail,
  nl.
go => true.

model(Data) =>
  Len = Data.len,

  BetaA = uniform(0,Len),
  BetaB = uniform(0,Len),
  Prob  = beta_dist(BetaA,BetaB),

  IsRightHanded = flip_n(Prob,Len),
  foreach(I in 1..Len)
    observe(IsRightHanded[I] == Data[I])
  end,

  if observed_ok then
    add("betaA",BetaA),
    add("betaB",BetaB),
    add("prob",Prob)
  end.