/* 

  Multinomial voting in Picat.

  Multinomial example voting
  From https://en.wikipedia.org/wiki/Multinomial_distribution#Example
  """
  Suppose that in a three-way election for a large country, candidate A received 20% of 
  the votes, candidate B received 30% of the votes, and candidate C received 50% of 
  the votes. 

  If six voters are selected randomly, what is the probability that there will be exactly 
  one supporter for candidate A, two supporters for candidate B and three supporters for 
  candidate C in the sample?

  Note: Since we’re assuming that the voting population is large, it is reasonable 
  and permissible to think of the probabilities as unchanging once a voter is selected 
  for the sample. Technically speaking this is sampling without replacement, so the 
  correct distribution is the multivariate hypergeometric distribution, but the 
  distributions converge as the population grows large.

  Pr = ... = 0.135 
   """

  Cf my Gamble model gamble_multinomial_voting.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.

/*
  Using categorical_dist

  Probabilities:
  [1,2,3]: 0.1335000000000000
  [1,1,4]: 0.1113000000000000
  [2,1,3]: 0.0895000000000000
  [2,2,2]: 0.0872000000000000
  [0,2,4]: 0.0859000000000000
  [1,3,2]: 0.0834000000000000
  [0,3,3]: 0.0642000000000000
  [0,1,5]: 0.0580000000000000
  [2,0,4]: 0.0389000000000000
  [1,0,5]: 0.0369000000000000
  [3,1,2]: 0.0359000000000000
  [0,4,2]: 0.0307000000000000
  [2,3,1]: 0.0300000000000000
  [1,4,1]: 0.0229000000000000
  [3,2,1]: 0.0215000000000000
  [3,0,3]: 0.0203000000000000
  [0,0,6]: 0.0140000000000000
  [4,1,1]: 0.0088000000000000
  [0,5,1]: 0.0065000000000000
  [2,4,0]: 0.0056000000000000
  [4,0,2]: 0.0045000000000000
  [3,3,0]: 0.0045000000000000
  [1,5,0]: 0.0026000000000000
  [4,2,0]: 0.0023000000000000
  [5,0,1]: 0.0006000000000000
  [0,6,0]: 0.0004000000000000
  [5,1,0]: 0.0001000000000000

  var : p
  Probabilities:
  false: 0.8665000000000000
  true: 0.1335000000000000

  Theoretical: 0.13500000000000


*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs
                           % show_percentiles,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  printf("Theoretical: %.14f\n",multinomial_dist_pdf(6,[20/100,30/100,50/100],[1,2,3])),
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  % Percentage of votes for A, B, and C
  Pct = [20/100,30/100,50/100],

  % There are 6 votes
  M = categorical_dist_n(Pct,[a,b,c],6),
  Os = count_occurrences_list([a,b,c],M),

  % Probability of 1 vote for A, 2 votes for B, and 3 votes for C
  P = check(Os == [1,2,3]),
  
  add("os",Os),
  add("p",P).


/*
  Using multinomial_dist for simulation

  var : m
  Probabilities:
  [1,2,3]: 0.1358000000000000
  [1,1,4]: 0.1133000000000000
  [2,1,3]: 0.0894000000000000
  [0,2,4]: 0.0842000000000000
  [1,3,2]: 0.0825000000000000
  [2,2,2]: 0.0767000000000000
  [0,3,3]: 0.0686000000000000
  [0,1,5]: 0.0604000000000000
  [1,0,5]: 0.0398000000000000
  [2,0,4]: 0.0372000000000000
  [3,1,2]: 0.0336000000000000
  [2,3,1]: 0.0334000000000000
  [0,4,2]: 0.0288000000000000
  [1,4,1]: 0.0260000000000000
  [3,2,1]: 0.0212000000000000
  [3,0,3]: 0.0201000000000000
  [0,0,6]: 0.0149000000000000
  [4,1,1]: 0.0070000000000000
  [0,5,1]: 0.0069000000000000
  [4,0,2]: 0.0049000000000000
  [2,4,0]: 0.0042000000000000
  [3,3,0]: 0.0040000000000000
  [4,2,0]: 0.0024000000000000
  [1,5,0]: 0.0021000000000000
  [0,6,0]: 0.0012000000000000
  [5,0,1]: 0.0011000000000000
  [5,1,0]: 0.0003000000000000

  var : p
  Probabilities:
  false: 0.8642000000000000
  true: 0.1358000000000000

  Theoretical: 0.13500000000000

*/
go2 ?=>
  reset_store,
  run_model(10_000,$model2,[show_probs
                           % show_percentiles,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  printf("Theoretical: %.14f\n",multinomial_dist_pdf(6,[20/100,30/100,50/100],[1,2,3])),
  % show_store_lengths,nl,
  % fail,
  nl.
go2 => true.
  
  
model2() =>
  % Percentage of votes for A, B, and C
  Pct = [20/100,30/100,50/100],

  % There are 6 votes
  M = multinomial_dist(6,Pct),

  % Probability of 1 vote for A, 2 votes for B, and 3 votes for C
  P = check(M == [1,2,3]),
  
  add("m",M),
  add("p",P).