/* 

  Three way election in Picat.

  From Statistics101 (Resampling Stats)
  File threewayElection.txt
  """
  A poll was taken before an election in which three politicians were
  competing. The sample size was 400. The result was that 10 percent of
  those polled favored politician A, 40 percent favored politician B,
  and the rest favored politician C. What are the 95 percent confidence
  intervals for the results for each politician?
  -> 
  polAConfidence: (7.249999999999999 13.0)
  polBConfidence: (35.25 44.75)
  polCConfidence: (45.25 54.75)
  """

  Two models below:
  - model: resampling
  - model2: multinomial


  Cf my Gamble model gamble_three_way_election.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.
% import ordset.

main => go.

/*
  Resampling

  var : a
  mean = 9.99393
  HPD intervals:
  HPD interval (0.95): 7.00000000000000..12.75000000000000

  var : b
  mean = 39.9988
  HPD intervals:
  HPD interval (0.95): 35.00000000000000..44.50000000000000

  var : c
  mean = 50.007275
  HPD intervals:
  HPD interval (0.95): 44.75000000000000..54.50000000000000

  var : a > 50
  mean = 0.0

  var : b > 50
  mean = 0.0

  var : c > 50
  mean = 0.4801

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[mean,
                           show_hpd_intervals,hpd_intervals=[0.95]
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  Population = rep(10,a) ++ rep(40,b) ++ rep(50,c),
  N = 400,
  Sample = resample(N,Population),

  A = 100 * count_occurrences(a,Sample) / N,
  B = 100 * count_occurrences(b,Sample) / N,
  C = 100 * count_occurrences(c,Sample) / N,

  add("a",A),
  add("b",B),
  add("c",C),
  add("a > 50",check(A > 50)),
  add("b > 50",check(B > 50)),
  add("c > 50",check(C > 50)).


/*
  Using multinomial

  var : a
  mean = 9.99675
  HPD intervals:
  HPD interval (0.95): 7.00000000000000..12.75000000000000

  var : b
  mean = 39.9969
  HPD intervals:
  HPD interval (0.95): 35.00000000000000..44.50000000000000

  var : c
  mean = 50.0063
  HPD intervals:
  HPD interval (0.95): 45.25000000000000..54.75000000000000

  var : a > 50
  mean = 0.0

  var : b > 50
  mean = 0.0

  var : c > 50
  mean = 0.4795


*/
go2 ?=>
  reset_store,
  run_model(10_000,$model2,[mean,
                           show_hpd_intervals,hpd_intervals=[0.95]
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go2 => true.
  
  
model2() =>
  N = 400,
  Sample = multinomial_dist(N,[10/100,40/100,50/100]),

  A = 100 * Sample[1] / N,
  B = 100 * Sample[2] / N,
  C = 100 * Sample[3] / N,

  add("a",A),
  add("b",B),
  add("c",C),
  add("a > 50",check(A > 50)),
  add("b > 50",check(B > 50)),
  add("c > 50",check(C > 50)).