/* 

  Rademacher dist in Picat.

  https://en.wikipedia.org/wiki/Rademacher_distribution
  """
  In probability theory and statistics, the Rademacher distribution (which is named after 
  Hans Rademacher) is a discrete probability distribution where a random variate X has a 50% 
  chance of being +1 and a 50% chance of being -1.
  """


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

/*

  var : x
  Probabilities:
  -1: 0.5080000000000000
  1: 0.4920000000000000
  mean = -0.016

  var : sum y
  Probabilities:
  0: 0.2424000000000000
  -2: 0.2131000000000000
  2: 0.1996000000000000
  -4: 0.1223000000000000
  4: 0.1158000000000000
  -6: 0.0449000000000000
  6: 0.0414000000000000
  -8: 0.0097000000000000
  8: 0.0090000000000000
  10: 0.0009000000000000
  -10: 0.0009000000000000
  mean = -0.0796

  var : last y
  Probabilities:
  1: 0.5024999999999999
  -1: 0.4975000000000000
  mean = 0.005

  var : a
  Probabilities:
  0: 0.2424000000000000
  -2: 0.2131000000000000
  2: 0.1996000000000000
  -4: 0.1223000000000000
  4: 0.1158000000000000
  -6: 0.0449000000000000
  6: 0.0414000000000000
  -8: 0.0097000000000000
  8: 0.0090000000000000
  10: 0.0009000000000000
  -10: 0.0009000000000000
  mean = -0.0796


  rand = [1,-1,1,-1,-1,-1,-1,1,-1,-1]
  pdf = 0.5
  cdf = 1
  quantile = 1


*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs,mean % ,
                           % show_percentiles,show_histogram,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  println(rand=rademacher_dist_n(10)),  
  println(pdf=rademacher_dist_pdf(1)),
  println(cdf=rademacher_dist_cdf(1)),
  println(quantile=rademacher_dist_quantile(1,0.99)),    
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.

model() =>
  X = rademacher_dist(),
  Y = rademacher_dist_n(10),
  SumY = Y.sum,
  LastY = Y.last,
  A = Y.cumsum.last,
  
  add("x",X),
  add("sum y",SumY),
  add("last y",LastY),
  add("a",A).