/* 

  Logit distribution in Picat.

  From https://en.wikipedia.org/wiki/Logit
  """
  In statistics, the logit function is the quantile function associated with the standard
  logistic distribution. It has many uses in data analysis and machine learning, 
  especially in data transformations.

  Mathematically, the logit is the inverse of the standard logistic function 
  sigma (x)=1/(1+e^(-x)), so the logit is defined as
    logit(p) = sigma^-1(p) = ln(p/(1-p))  p in 0..1
   """

  Note: According to Mathematica LogisticDistribution, the quantile of 
  LogisticDistribution[Mu,Beta,x] is
    Mu - Beta*Log[-1 + 1/x]


  Cf my Gamble model gamble_logit_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 : d
Probabilities (truncated):
9.288769606217127: 0.0001000000000000
9.198537700995011: 0.0001000000000000
8.904182202923693: 0.0001000000000000
8.262305181805095: 0.0001000000000000
.........
-8.015821951305064: 0.0001000000000000
-8.11986515367056: 0.0001000000000000
-8.398195754223158: 0.0001000000000000
-9.167655205552682: 0.0001000000000000
mean = 0.00837819
Percentiles:
(0.001 -6.3799033891465475)
(0.01 -4.5299706351433597)
(0.025 -3.6006173923380098)
(0.05 -2.8987406249798466)
(0.1 -2.1647686072779546)
(0.25 -1.0740296008823442)
(0.5 -0.00097435647128925405)
(0.75 1.1042650759272994)
(0.84 1.6522016510397421)
(0.9 2.2157671814910005)
(0.95 2.9120094058036279)
(0.975 3.6279359969930658)
(0.99 4.6477836105512296)
(0.999 6.8030195906804725)
(0.9999 9.1985467241854693)
(0.99999 9.2797473180140759)
HPD intervals:
HPD interval (0.5): -1.01472422579231..1.14949890125836
HPD interval (0.84): -2.30240791254364..2.51975624067660
HPD interval (0.95): -3.47334256229409..3.72780852375865
Histogram (total 10000)
 -9.168:     1                                                              (0.000 / 0.000)
 -8.706:     0                                                              (0.000 / 0.000)
 -8.245:     3                                                              (0.000 / 0.000)
 -7.783:     1                                                              (0.000 / 0.001)
 -7.322:     1                                                              (0.000 / 0.001)
 -6.861:     3                                                              (0.000 / 0.001)
 -6.399:     4                                                              (0.000 / 0.001)
 -5.938:    13 #                                                            (0.001 / 0.003)
 -5.476:    15 #                                                            (0.002 / 0.004)
 -5.015:    39 ##                                                           (0.004 / 0.008)
 -4.554:    49 ###                                                          (0.005 / 0.013)
 -4.092:    63 ###                                                          (0.006 / 0.019)
 -3.631:   118 ######                                                       (0.012 / 0.031)
 -3.169:   173 #########                                                    (0.017 / 0.048)
 -2.708:   268 ##############                                               (0.027 / 0.075)
 -2.246:   404 #####################                                        (0.040 / 0.115)
 -1.785:   565 ##############################                               (0.057 / 0.172)
 -1.324:   750 #######################################                      (0.075 / 0.247)
 -0.862:   979 ###################################################          (0.098 / 0.345)
 -0.401:  1128 ###########################################################  (0.113 / 0.458)
  0.061:  1146 ############################################################ (0.115 / 0.572)
  0.522:  1071 ########################################################     (0.107 / 0.679)
  0.983:   922 ################################################             (0.092 / 0.772)
  1.445:   707 #####################################                        (0.071 / 0.842)
  1.906:   519 ###########################                                  (0.052 / 0.894)
  2.368:   390 ####################                                         (0.039 / 0.933)
  2.829:   226 ############                                                 (0.023 / 0.956)
  3.290:   158 ########                                                     (0.016 / 0.972)
  3.752:   108 ######                                                       (0.011 / 0.982)
  4.213:    61 ###                                                          (0.006 / 0.989)
  4.675:    42 ##                                                           (0.004 / 0.993)
  5.136:    26 #                                                            (0.003 / 0.995)
  5.597:    17 #                                                            (0.002 / 0.997)
  6.059:    14 #                                                            (0.001 / 0.998)
  6.520:     5                                                              (0.001 / 0.999)
  6.982:     4                                                              (0.000 / 0.999)
  7.443:     2                                                              (0.000 / 1.000)
  7.905:     1                                                              (0.000 / 1.000)
  8.366:     1                                                              (0.000 / 1.000)
  8.827:     3                                                              (0.000 / 1.000)


*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean,
                           show_percentiles,
                           show_hpd_intervals,hpd_intervals=[0.5,0.84,0.95],
                           show_histogram
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  

logit(P) = log(P/(1-P)).
logit_dist() = logit(U) =>
  U = uniform(0,1).

model() =>
  
  D  = logit_dist(),
  add("d",D).