/* 

  Pig food in Picat.

  From "Resampling Stats Illustrations" (http://www.statistics101.net/PeterBruce_05-illus.pdf)
  Page 29
  """
  Pig weight gains — reliability of the estimate
  (program “pigfood”)
  (This bootstrap example is from Basic Research Methods in Social
  Science, Julian L. Simon, 1969.)
  An agricultural lab decides to experiment with a new pig ration —
  ration A — on twelve pigs. After 4 weeks, the pigs experience an
  average gain of 508 ounces. The weight gains of the individual
  pigs are as follows: 496, 544, 464, 416, 512, 560, 608, 544, 480, 466,
  512, 496.
  In presenting these results to major agricultural feed distributors,
  the lab wants to report not just the average estimated weight gain
  (as represented by the sample average), but also the possible range
  of sampling error.
  How can we determine the extent to which one sample differs
  from another? (The reliability of our estimated mean weight gain.)
  If we had more time and money, we could try the ration on addi-
  tional groups of 12 pigs, and see how the mean weight gain
  differed from group to group.
  -> 
  interval = 480 to 537 
  """

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

/*
  Here we use some different HPD-intervals

  var : samples mean
  Probabilities (truncated):
  508.166666666666686: 0.0166000000000000
  509.5: 0.0165000000000000
  505.5: 0.0152000000000000
  502.833333333333314: 0.0148000000000000
  .........
  463.0: 0.0001000000000000
  461.666666666666686: 0.0001000000000000
  461.5: 0.0001000000000000
  461.333333333333314: 0.0001000000000000
  mean = 508.089
  HPD intervals:
  HPD interval (0.75): 492.16666666666669..524.16666666666663
  HPD interval (0.84): 488.16666666666669..526.83333333333337
  HPD interval (0.9): 485.50000000000000..530.83333333333337
  HPD interval (0.95): 480.16666666666669..534.66666666666663
  HPD interval (0.975): 477.66666666666669..540.16666666666663
  HPD interval (0.99): 472.16666666666669..544.00000000000000

*/
go ?=>
  Weights = [496,544,464,416,512,560,608,544,480,466,512,496],
  reset_store,
  run_model(10_000,$model(Weights),[show_probs_trunc,mean,
                           % show_percentiles,
                           show_hpd_intervals,hpd_intervals=[0.75,0.84,0.9,0.95,0.975, 0.99]
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model(Weights) =>
  Samples = resample(Weights),
  SamplesMean = Samples.mean,
  
  add("samples mean",SamplesMean).