#| 

  Min stable distribution in Racket/Gamble 

  From Mathematica's MinStableDistribution

  This is also called Generalized Extreme Value distribution 
  (for minimum values).

  This program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Racket page: http://www.hakank.org/racket/

|#

#lang gamble

; (require gamble/viz)
(require racket)
(require "gamble_utils.rkt")
(require "gamble_distributions.rkt")


#|
  variable : d
  2.37888723606864: 0.0009999999999999994
  2.825076522317268: 0.0009999999999999994
  2.4921439654515227: 0.0009999999999999994
  1.9320985264197237: 0.0009999999999999994
  2.011445995191115: 0.0009999999999999994
  ...
  1.5019607984432175: 0.0009999999999999994
  2.496320875502025: 0.0009999999999999994
  0.8917664696369227: 0.0009999999999999994
  -0.29066928446632856: 0.0009999999999999994
  -3.530076377127445: 0.0009999999999999994
  mean: -1.4867262687602532

  variable : p
  #t: 0.6480000000000005
  #f: 0.35200000000000026
  mean: 0.6480000000000005

|#
(define (model)
  (; enumerate
   ; rejection-sampler
   importance-sampler
   ; mh-sampler

   (define mu 2)
   (define sigma 1)
   (define xi 0.91)

   (define d (min_stable_dist mu sigma xi))

   (define p (<= d mu))
   
   (list d
         p)

   )
)

(show-marginals (model)
                (list  "d"
                       "p"
                       )
                #:num-samples 1000
                #:truncate-output 5
                ; #:skip-marginals? #t
                ; #:show-stats? #t
                ; #:credible-interval 0.84
                ; #:hpd-interval (list 0.84)
                ; #:show-histogram? #t
                ; #:show-percentiles? #t
                ; #:burn 0
                ; #:thin 0
                )

#|
  Recovering parameters

  For 100 values of min (normal 100 15)

  data: (49.56433300638627 45.42849077209317 51.12661125815772 50.09879927148479 45.50982545788071 50.12364736831026 45.978539150875996 49.65552749984283 52.67266622380643 48.56863948649306 54.231166133642915 59.065258163946545 45.6658753081055 31.870492901612224 50.53267786344057 52.18714001944357 55.27130731717705 56.28023670315283 46.85955601292287 53.9609220940334 58.2374921389876 49.596799437032715 54.85320532843294 49.2321506009014 45.762462316057146 57.59501633584913 45.03006869798502 52.68256202513432 51.23253490640521 56.98248810563077 53.61508467411477 46.421500086857044 52.745147986087986 48.908633087978046 50.71634927583853 56.79605748003404 53.02904735697946 50.502628511092546 51.98423720882249 48.808798570154885 51.4258793859458 46.58559344077984 44.07384898797801 53.82228545922746 41.91380122134164 52.81367781779215 49.91672485936912 45.89670531267542 49.225072397970436 52.04107249106037 53.95381609587204 57.34502197561476 60.529307085528316 56.55769045913069 49.57077310532649 46.53700697366534 45.902226178340456 51.603037375657124 43.37475734284823 49.56486229202943 58.29593276733964 44.29571289319839 45.727321817759744 46.51637108053079 50.56938904476137 56.0741194820407 47.821823077124094 50.91493687042285 56.4832165418988 51.80865295374938 61.12521242458355 53.97978534731618 49.72069768389238 55.6680554836353 33.148481544576754 56.833468734329095 50.89487469832344 53.277366930146655 61.669533741147426 51.54461809832464 48.841573439761255 57.89393195499958 44.06855830745549 53.754390014959476 56.45422928049421 53.31680731650925 55.33044567596348 58.80820670138492 55.06711052173423 48.76597988613222 50.302001654447366 59.21120594301467 58.96156160598758 49.550544006863305 55.482698939983 47.92797168785546 51.631633976910166 56.11236483463039 57.610370764861074 47.31249472825772)
(min 31.870492901612224 mean 51.34810788858642 max 61.669533741147426 stddev 5.207059104612523 variance 27.11346451892817)

  31.87:  1 #### (0.01  / 0    )
  34.85:  1 #### (0.01  / 0.01 )
  37.83:  0  (0     / 0.02 )
  40.81:  0  (0     / 0.02 )
  43.79:  2 ######## (0.02  / 0.02 )
  46.77: 16 ########################################################### (0.16  / 0.04 )
  49.75: 18 ################################################################## (0.18  / 0.2  )
  52.73: 22 ################################################################################ (0.22  / 0.38 )
  55.71: 18 ################################################################## (0.18  / 0.6  )
  58.69: 15 ####################################################### (0.15  / 0.78 )

  variable : mu
  mean: 51.42012329201724
  HPD interval (0.001): 48.416897083675714..48.585113000576214
  HPD interval (0.84): 50.06642295473801..52.939660286458995
  HPD interval (0.999): 48.24638012528553..54.526635399588855

  variable : sigma
  mean: 0.8087017741601884
  HPD interval (0.001): 0.0038495056332780915..0.0038495056332780915
  HPD interval (0.84): 0.0038495056332780915..1.571139638963399
  HPD interval (0.999): 0.0038495056332780915..3.6777341936176446

  variable : xi
  mean: -0.881649123866001
  HPD interval (0.001): -2.8893365922807743..-2.8893365922807743
  HPD interval (0.84): -2.136688320066587..0.20465945535552843
  HPD interval (0.999): -2.8893365922807743..1.8435837362113916

  variable : post
  mean: 51.24724480988478
  HPD interval (0.001): 43.99200050866426..44.19101623721647
  HPD interval (0.84): 49.32439037898759..52.957057673835195
  HPD interval (0.999): 43.99200050866426..62.72244107101271

  variable : p
  mean: 8.964342171394439e-11


|#
(define data (for/list ([i 100])
               (min-list (for/list ([i 1000])
                           (normal 100 15)
                           ; (binomial 100 0.5)
                           ; (exponential 10)
                           ; (gamma 10 10)
                           ; (poisson 10)
                           ; (sample (t-dist 2 100 15)) ; The parameters are in different order from Mathematica's StudentTDistribution. This is slow!
                           ))))

(show "data" data)
(show2 "min" (min-list data) "mean" (avg data) "max" (max-list data) "stddev" (stddev data) "variance" (variance data))
(newline)
(show-histogram data)
(newline)
(define (model2)
  (; enumerate
   ; rejection-sampler
   importance-sampler
   ; mh-sampler

   ; (define mu (normal (avg data) 10))
   (define mu (normal (avg data) (stddev data)))
   ; (define mu (uniform (min-list data) (max-list data)))
   (define sigma (uniform 0 20))
   (define xi (uniform -3 3))

   (for ([i (length data)])
     (observe-sample (normal-dist (min_stable_dist mu sigma xi) 10) (list-ref data i)))
   
   (define post (min_stable_dist mu sigma xi))

   (define p (<= post (min-list data)))
   
   (list mu
         sigma
         xi
         post
         p
         )

   )
)

(show-marginals (model2)
                (list  "mu"
                       "sigma"
                       "xi"
                       "post"
                       "p"
                       )
                #:num-samples 10000
                #:truncate-output 5
                #:skip-marginals? #t
                ; #:show-stats? #t
                ; #:credible-interval 0.84
                #:hpd-interval (list 0.001 0.84 0.999)
                ; #:show-histogram? #t
                ; #:show-percentiles? #t
                ; #:burn 0
                ; #:thin 0
                )