/* 

  PPL - Geometric coin in Picat.

  From Mathematica (GeometricDistribution) 
  """
  A coin-tossing experiment consists of tossing a fair coin repeatedly 
  until a tail results. Simulate the process:
    RandomVariate(GeometricDistribution(1/2), 10)
    -> 
     (0,0,1,1,2,0,0,2,0,1)

  Compute the probability that at least 4 coin tosses will be necessary:
    Probability(x >= 4, x -> GeometricDistribution(1/2))
    -> 1/16 (0.0625)
    
  Compute the expected number of coin tosses:
    Mean(GeometricDistribution(1/2)) 
    -> 1

  """

  Picat> X=geometric_dist_n(1/2,10)
  X = [0,3,2,0,0,1,3,5,1,3]

  Picat> X=1-geometric_dist_cdf(1/2,3)
  X = 0.0625

  Picat> X=geometric_dist_mean(1/2)   
  X = 1.0


  This is a port of my Racket/Gamble model gamble_geometric_coin.pi

  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.


/*
Probabilities:
0: 0.4981000000000000
1: 0.2545000000000000
2: 0.1243000000000000
3: 0.0608000000000000
4: 0.0307000000000000
5: 0.0170000000000000
6: 0.0075000000000000
7: 0.0033000000000000
8: 0.0018000000000000
9: 0.0011000000000000
10: 0.0006000000000000
11: 0.0002000000000000
12: 0.0001000000000000
mean = 0.9951
Percentiles:
(0.001 0)
(0.01 0)
(0.025 0)
(0.05 0)
(0.1 0)
(0.25 0)
(0.5 1)
(0.75 1)
(0.84 2)
(0.9 3)
(0.95 4)
(0.975 5)
(0.99 6)
(0.999 9)
(0.9999 11.000099999999293)
(0.99999 11.900010000001203)
HPD intervals:
HPD interval (0.5): 0.00000000000000..1.00000000000000
HPD interval (0.84): 0.00000000000000..2.00000000000000
HPD interval (0.9): 0.00000000000000..3.00000000000000
HPD interval (0.95): 0.00000000000000..4.00000000000000
HPD interval (0.99): 0.00000000000000..6.00000000000000
HPD interval (0.999): 0.00000000000000..9.00000000000000
HPD interval (0.9999): 0.00000000000000..11.00000000000000
HPD interval (0.99999): 0.00000000000000..12.00000000000000
Histogram:
      0:  4981 ############################################################ (0.498 / 0.498)
      1:  2545 ###############################                              (0.255 / 0.753)
      2:  1243 ###############                                              (0.124 / 0.877)
      3:   608 #######                                                      (0.061 / 0.938)
      4:   307 ####                                                         (0.031 / 0.968)
      5:   170 ##                                                           (0.017 / 0.985)
      6:    75 #                                                            (0.007 / 0.993)
      7:    33                                                              (0.003 / 0.996)
      8:    18                                                              (0.002 / 0.998)
      9:    11                                                              (0.001 / 0.999)
     10:     6                                                              (0.001 / 1.000)
     11:     2                                                              (0.000 / 1.000)
     12:     1                                                              (0.000 / 1.000)

Var p_at_least_4:
Probabilities:
0: 0.9377000000000000
1: 0.0623000000000000
mean = 0.0623
Percentiles:
(0.001 0)
(0.01 0)
(0.025 0)
(0.05 0)
(0.1 0)
(0.25 0)
(0.5 0)
(0.75 0)
(0.84 0)
(0.9 0)
(0.95 1)
(0.975 1)
(0.99 1)
(0.999 1)
(0.9999 1)
(0.99999 1)
HPD intervals:
HPD interval (0.5): 0.00000000000000..0.00000000000000
HPD interval (0.84): 0.00000000000000..0.00000000000000
HPD interval (0.9): 0.00000000000000..0.00000000000000
HPD interval (0.95): 0.00000000000000..1.00000000000000
HPD interval (0.99): 0.00000000000000..1.00000000000000
HPD interval (0.999): 0.00000000000000..1.00000000000000
HPD interval (0.9999): 0.00000000000000..1.00000000000000
HPD interval (0.99999): 0.00000000000000..1.00000000000000
Histogram:
      0:  9377 ############################################################ (0.938 / 0.938)
      1:   623 ####                                                         (0.062 / 1.000)


using_cdf = 0.0625 (1 / 16)

samples = [1,0,0,0,1,1,2,3,0,1,0,2,0,0,1,0,2,1,1,1]

*/

go ?=>
  reset_store(),
  run_model(10000,geometric_coin,),
  ExactProb = 1-geometric_dist_cdf(1/2,3),
  println(using_cdf=ExactProb.and_rat),
  nl,
  println(samples=[geometric_dist(1/2) : _ in 1..20]),
  nl.
go => true.

geometric_coin() =>
  V = geometric_dist(1/2),
  P_at_least_4 = cond(V >= 4,1,0),
  add_all([["p_at_least_4",P_at_least_4],
           ["V",V]]).