/* 

  Profits in Picat.

  From "Resampling Stats Illustrations" (http://www.statistics101.net/PeterBruce_05-illus.pdf)
  Page 33
  """
  A confidence interval for a median
  (program “profits”)
  Sometimes, especially with highly skewed data such as incomes,
  the median is preferred to the mean as a measure of the distribu-
  tion center. Constructing a confidence interval for the median is
  easy with RESAMPLING STATS.
  Say you need to come up with a quick estimate of the profits of a
  typical American Fortune 1000 business, and the extent to which
  that estimate might be in error. You draw a random sample of 15
  firms, finding their profits (in $ million) to be: 1315, 288, 155, 37,
  99, 40, 170, 66, 500, 419, 125, -90, -63, 29, 966. We use
  RESAMPLING STATS to calculate the median profit, and construct
  a bootstrap confidence interval
  -> 
  median = 125
  interval = 40 288
  """

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

/*

  profits = [1315,288,155,37,99,40,170,66,500,419,125,-90,-63,29,966]
  [min = -90,mean = 270.4,median = 125,max = 1315,variance = 144849.0,stdev = 380.59]

  var : sample median
  mean = 135.082
  HPD intervals:
  HPD interval (0.9): 37.00000000000000..170.00000000000000
  HPD interval (0.95): 40.00000000000000..288.00000000000000
  HPD interval (0.99): 37.00000000000000..419.00000000000000
  Percentiles:
  (0.001 29)
  (0.01 37)
  (0.025 40)
  (0.05 40)
  (0.1 66)
  (0.25 99)
  (0.5 125)
  (0.75 155)
  (0.84 170)
  (0.9 170)
  (0.95 288)
  (0.975 419)
  (0.99 419)
  (0.999 500)
  (0.9999 966)
  (0.99999 966)

  var : sample mean
  mean = 270.863
  HPD intervals:
  HPD interval (0.9): 104.73333333333333..425.13333333333333
  HPD interval (0.95): 95.33333333333333..470.66666666666669
  HPD interval (0.99): 55.33333333333334..535.46666666666670
  Percentiles:
  (0.001 37.531999999999996)
  (0.01 77.064666666666668)
  (0.025 99.594999999999999)
  (0.05 120.93333333333334)
  (0.1 147.59999999999999)
  (0.25 201.26666666666668)
  (0.5 264.36666666666667)
  (0.75 333.39999999999998)
  (0.84 370.60000000000002)
  (0.9 402.94000000000005)
  (0.95 448.13333333333333)
  (0.975 479.87499999999989)
  (0.99 521.13866666666672)
  (0.999 603.53726666666739)
  (0.9999 656.28221333322335)
  (0.99999 796.18822133352035)

  var : sample min
  mean = -69.038
  HPD intervals:
  HPD interval (0.9): -90.00000000000000..29.00000000000000
  HPD interval (0.95): -90.00000000000000..29.00000000000000
  HPD interval (0.99): -90.00000000000000..40.00000000000000
  Percentiles:
  (0.001 -90)
  (0.01 -90)
  (0.025 -90)
  (0.05 -90)
  (0.1 -90)
  (0.25 -90)
  (0.5 -90)
  (0.75 -63)
  (0.84 -63)
  (0.9 29)
  (0.95 29)
  (0.975 37)
  (0.99 40)
  (0.999 66)
  (0.9999 99.005599999960396)
  (0.99999 149.40056000006734)

  var : sample max
  mean = 1131.59
  HPD intervals:
  HPD interval (0.9): 500.00000000000000..1315.00000000000000
  HPD interval (0.95): 500.00000000000000..1315.00000000000000
  HPD interval (0.99): 419.00000000000000..1315.00000000000000
  Percentiles:
  (0.001 170)
  (0.01 419)
  (0.025 419)
  (0.05 500)
  (0.1 500)
  (0.25 966)
  (0.5 1315)
  (0.75 1315)
  (0.84 1315)
  (0.9 1315)
  (0.95 1315)
  (0.975 1315)
  (0.99 1315)
  (0.999 1315)
  (0.9999 1315)
  (0.99999 1315)



*/
go ?=>
  Profits = [1315,288,155,37,99,40,170,66,500,419,125,-90,-63,29,966],
  println(profits=Profits),
  show_simple_stats(Profits),
  reset_store,
  run_model(10_000,$model(Profits),[% show_probs_trunc,
                           mean,
                           show_hpd_intervals,hpd_intervals=[0.9,0.95,0.99],
                           show_percentiles                           
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model(Profits) =>
  Sample = resample(Profits),
  SampleMedian = Sample.median,
  SampleMean   = Sample.mean,
  SampleMin    = Sample.min,
  SampleMax    = Sample.max,
  
  
  add("sample median",SampleMedian),
  add("sample mean",SampleMean),
  add("sample min",SampleMin),
  add("sample max",SampleMax).