/* 

  Kumaraswamy distribution in Picat.

  From Handbook on probability distributions
  """
  Since the quantile function is explicit
     F^-1(u) = (1 - (1 - u)^(1/b))^(1/a)
  an inversion function method F^-1(u) with u uniformly distributed is easily computable.
  """

  Picat> X=kumaraswamy_dist_mean(5,2)
  X = 0.757575757575757


  Cf my Gamble model gamble_kumaraswamy_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 : g
Probabilities (truncated):
0.999754217446404: 0.0001000000000000
0.998635732458591: 0.0001000000000000
0.996333556885938: 0.0001000000000000
0.996259090745204: 0.0001000000000000
.........
0.193213083730812: 0.0001000000000000
0.184008116975508: 0.0001000000000000
0.12473335384345: 0.0001000000000000
0.110999359839025: 0.0001000000000000
mean = 0.757562
HPD intervals:
HPD interval (0.94): 0.49779165545513..0.99529796493627
Histogram (total 10000)
  0.111:     1                                                              (0.000 / 0.000)
  0.133:     1                                                              (0.000 / 0.000)
  0.155:     0                                                              (0.000 / 0.000)
  0.178:     1                                                              (0.000 / 0.000)
  0.200:     2                                                              (0.000 / 0.001)
  0.222:     6 #                                                            (0.001 / 0.001)
  0.244:     9 #                                                            (0.001 / 0.002)
  0.267:     8 #                                                            (0.001 / 0.003)
  0.289:    28 ##                                                           (0.003 / 0.006)
  0.311:    19 ##                                                           (0.002 / 0.007)
  0.333:    28 ##                                                           (0.003 / 0.010)
  0.355:    37 ###                                                          (0.004 / 0.014)
  0.378:    44 ####                                                         (0.004 / 0.018)
  0.400:    63 ######                                                       (0.006 / 0.025)
  0.422:    58 #####                                                        (0.006 / 0.030)
  0.444:    85 ########                                                     (0.009 / 0.039)
  0.467:   121 ###########                                                  (0.012 / 0.051)
  0.489:    97 #########                                                    (0.010 / 0.061)
  0.511:   160 ##############                                               (0.016 / 0.077)
  0.533:   162 ##############                                               (0.016 / 0.093)
  0.555:   193 #################                                            (0.019 / 0.112)
  0.578:   235 #####################                                        (0.024 / 0.136)
  0.600:   255 #######################                                      (0.025 / 0.161)
  0.622:   299 ###########################                                  (0.030 / 0.191)
  0.644:   299 ###########################                                  (0.030 / 0.221)
  0.666:   407 ####################################                         (0.041 / 0.262)
  0.689:   433 ######################################                       (0.043 / 0.305)
  0.711:   450 ########################################                     (0.045 / 0.350)
  0.733:   539 ################################################             (0.054 / 0.404)
  0.755:   569 ###################################################          (0.057 / 0.461)
  0.778:   582 ####################################################         (0.058 / 0.519)
  0.800:   633 ########################################################     (0.063 / 0.582)
  0.822:   611 ######################################################       (0.061 / 0.644)
  0.844:   622 #######################################################      (0.062 / 0.706)
  0.866:   675 ############################################################ (0.068 / 0.773)
  0.889:   588 ####################################################         (0.059 / 0.832)
  0.911:   559 ##################################################           (0.056 / 0.888)
  0.933:   511 #############################################                (0.051 / 0.939)
  0.955:   374 #################################                            (0.037 / 0.976)
  0.978:   236 #####################                                        (0.024 / 1.000)


*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean,
                           % show_percentiles,
                           show_hpd_intervals,hpd_intervals=[0.94],
                           show_histogram
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  A = 5,
  B = 2,
  G = kumaraswamy_dist(A,B),
  
  add("g",G).  


/*
  Parameter recovery of kumaraswamy_dist(5,2)

  * Using (default) statistics mean and stdev

    Num accepted samples: 1000  Total samples: 105191  (0.010%)

    var : a
    Probabilities (truncated):
    9.940006125145596: 0.0010000000000000
    9.04091692340603: 0.0010000000000000
    8.7894361964564: 0.0010000000000000
    8.656433417506717: 0.0010000000000000
    .........
    3.871752400352504: 0.0010000000000000
    3.850304672303751: 0.0010000000000000
    3.79238008528593: 0.0010000000000000
    3.641090080156498: 0.0010000000000000
    mean = 6.11389

    var : b
    Probabilities (truncated):
    5.425299594022939: 0.0010000000000000
    4.583872255419322: 0.0010000000000000
    4.37933829419284: 0.0010000000000000
    4.281956851543838: 0.0010000000000000
    .........
    1.195042822591515: 0.0010000000000000
    1.194316819689384: 0.0010000000000000
    1.152894935865372: 0.0010000000000000
    0.946823838035028: 0.0010000000000000
    mean = 2.29526


  * Using statistics mean,stdev,kurtosis,skewness

    Num accepted samples: 1000  Total samples: 26181  (0.038%)
    var : a
    Probabilities (truncated):
    9.992217672589337: 0.0010000000000000
    9.990846903897285: 0.0010000000000000
    9.986955900684444: 0.0010000000000000
    9.986861712526: 0.0010000000000000
    .........
    2.418002154360526: 0.0010000000000000
    2.41038646002318: 0.0010000000000000
    2.149531729523806: 0.0010000000000000
    2.117233471659586: 0.0010000000000000
    mean = 7.18925

    var : b
    Probabilities (truncated):
    9.977769412839677: 0.0010000000000000
    9.939959763749483: 0.0010000000000000
    9.917992133068848: 0.0010000000000000
    9.907000304733868: 0.0010000000000000
    .........
    0.755586544785456: 0.0010000000000000
    0.722247518702525: 0.0010000000000000
    0.616867372201228: 0.0010000000000000
    0.508205000743365: 0.0010000000000000
    mean = 4.27586

*/
go2 ?=>
  Data = kumaraswamy_dist_n(5,2,100),
  println(data=Data),
  reset_store,
  run_model(10_000,$model2(Data),[show_probs_trunc,mean,
                           % show_percentiles,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % show_histogram
                           min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  show_store_lengths,nl,
  % fail,
  nl.
go2 => true.
  
  
model2(Data) =>
  A = uniform(0.01,10),
  B = uniform(0.01,10),  
  X = kumaraswamy_dist_n(A,B,Data.len),

  % observe_abc(Data,X,1/10),
  observe_abc(Data,X,1,[mean,stdev,kurtosis,skewness]),  

  if observed_ok then
    add("a",A), 
    add("b",B),   
  end.