/* 

  Earthquakes (Pareto dist) in Picat.

  From Mathematica ParetoDistribution
  """
  Consider earthquake magnitudes recorded in the US from 1935 to 1989:

  The integer parts of the magnitudes recorded on a Richter scale can be modeled with a 
  ParetoDistribution:
    edist = EstimatedDistribution[magnitudes, ParetoDistribution(k, Alpha, Gamma, Mu)
    -> 
    ParetoDistribution[46.1666, 4374., 0.326152, 2.]

  Find the probability of an earthquake with magnitude at least 6 on the Richter scale:
    NProbability[x >= 6, x -> edist]
    -> 0.0889446

  Find the average magnitude:
  Mean(edist)
  -> 4.68

  """

  For the 2 parameter version ParetoDistribution[k, alpha]:
  """
  edist = EstimatedDistribution[magnitudes,  ParetoDistribution[k, alpha]]
  -> ParetoDistribution[2., 1.18656]
  
  NProbability[x >= 6, x -> edist]
  """

  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.

/*

  For pareto1(K,Alpha):

[len = 2841,min = 2,mean = 4.7251,median = 5,max = 8,variance = 0.761105,stdev = 0.872414]
Histogram (total 2841)
      2:     5                                                              (0.002 / 0.002)
      3:    93 #####                                                        (0.033 / 0.034)
      4:  1208 ############################################################ (0.425 / 0.460)
      5:   965 ################################################             (0.340 / 0.799)
      6:   516 ##########################                                   (0.182 / 0.981)
      7:    50 ##                                                           (0.018 / 0.999)
      8:     4                                                              (0.001 / 1.000)
HPD intervals:
HPD interval (0.84): 4.00000000000000..6.00000000000000
HPD interval (0.94): 4.00000000000000..6.00000000000000
HPD interval (0.99): 3.00000000000000..7.00000000000000
methods = [mean,stdev,q25,q75]
explain = Data shows mild skewness or moderate tails. Adding quantiles for shape.
...
Num accepted samples: 100  Total samples: 1261  (0.079%)
var : k
Probabilities (truncated):
4.727095469239678: 0.0100000000000000
4.66776922483359: 0.0100000000000000
4.658383160437635: 0.0100000000000000
4.635091558487662: 0.0100000000000000
.........
3.243306063135763: 0.0100000000000000
3.236127805447266: 0.0100000000000000
3.179595839455536: 0.0100000000000000
3.080265425741796: 0.0100000000000000
mean = 3.98636
HPD intervals:
HPD interval (0.84): 3.52804808524812..4.62629930960308
HPD interval (0.94): 3.30254077045365..4.66776922483359
HPD interval (0.99): 3.17959583945554..4.72709546923968

var : alpha
Probabilities (truncated):
9.927869424516274: 0.0100000000000000
9.838867119810947: 0.0100000000000000
9.742890828588461: 0.0100000000000000
9.736620284494302: 0.0100000000000000
.........
4.163100151281385: 0.0100000000000000
3.940967501346473: 0.0100000000000000
3.933427106045852: 0.0100000000000000
3.809256382058029: 0.0100000000000000
mean = 7.07269
HPD intervals:
HPD interval (0.84): 4.93780948516811..9.60597900953422
HPD interval (0.94): 4.60102222967009..9.92786942451627
HPD interval (0.99): 3.93342710604585..9.92786942451627

var : post
Probabilities (truncated):
8.668419026784628: 0.0100000000000000
6.936818889746316: 0.0100000000000000
6.855458739199094: 0.0100000000000000
6.531924087241652: 0.0100000000000000
.........
3.480837847209461: 0.0100000000000000
3.453002742863971: 0.0100000000000000
3.416656379206965: 0.0100000000000000
3.186050953813038: 0.0100000000000000
mean = 4.6242
HPD intervals:
HPD interval (0.84): 3.57599107056948..5.47111562197377
HPD interval (0.94): 3.45300274286397..6.43912826615226
HPD interval (0.99): 3.18605095381304..6.93681888974632

var : p
Probabilities:
false: 0.9100000000000000
true: 0.0900000000000000
mean = 0.09

Post:

[len = 100,min = 3.18605,mean = 4.6242,median = 4.37685,max = 8.66842,variance = 0.79887,stdev = 0.893795]
Histogram (total 100)
  3.186:     1 ######                                                       (0.010 / 0.010)
  3.323:     0                                                              (0.000 / 0.010)
  3.460:     3 ##################                                           (0.030 / 0.040)
  3.597:     5 ##############################                               (0.050 / 0.090)
  3.734:     4 ########################                                     (0.040 / 0.130)
  3.871:     5 ##############################                               (0.050 / 0.180)
  4.008:    10 ############################################################ (0.100 / 0.280)
  4.145:    10 ############################################################ (0.100 / 0.380)
  4.283:    10 ############################################################ (0.100 / 0.480)
  4.420:     6 ####################################                         (0.060 / 0.540)
  4.557:     8 ################################################             (0.080 / 0.620)
  4.694:     5 ##############################                               (0.050 / 0.670)
  4.831:     6 ####################################                         (0.060 / 0.730)
  4.968:     2 ############                                                 (0.020 / 0.750)
  5.105:     2 ############                                                 (0.020 / 0.770)
  5.242:     5 ##############################                               (0.050 / 0.820)
  5.379:     4 ########################                                     (0.040 / 0.860)
  5.516:     2 ############                                                 (0.020 / 0.880)
  5.653:     1 ######                                                       (0.010 / 0.890)
  5.790:     1 ######                                                       (0.010 / 0.900)
  5.927:     1 ######                                                       (0.010 / 0.910)
  6.064:     1 ######                                                       (0.010 / 0.920)
  6.201:     0                                                              (0.000 / 0.920)
  6.338:     2 ############                                                 (0.020 / 0.940)
  6.475:     3 ##################                                           (0.030 / 0.970)
  6.613:     0                                                              (0.000 / 0.970)
  6.750:     0                                                              (0.000 / 0.970)
  6.887:     2 ############                                                 (0.020 / 0.990)
  7.024:     0                                                              (0.000 / 0.990)
  7.161:     0                                                              (0.000 / 0.990)
  7.298:     0                                                              (0.000 / 0.990)
  7.435:     0                                                              (0.000 / 0.990)
  7.572:     0                                                              (0.000 / 0.990)
  7.709:     0                                                              (0.000 / 0.990)
  7.846:     0                                                              (0.000 / 0.990)
  7.983:     0                                                              (0.000 / 0.990)
  8.120:     0                                                              (0.000 / 0.990)
  8.257:     0                                                              (0.000 / 0.990)
  8.394:     0                                                              (0.000 / 0.990)
  8.531:     1 ######                                                       (0.010 / 1.000)
HPD intervals:
HPD interval (0.84): 3.57599107056948..5.47111562197377
HPD interval (0.94): 3.45300274286397..6.43912826615226
HPD interval (0.99): 3.18605095381304..6.93681888974632


*/
go ?=>
  Magnitudes = [
5,5,4,6,7,4,4,5,4,5,4,4,6,4,5,5,5,6,5,4,4,6,6,
4,5,3,5,5,4,4,5,4,4,4,4,4,6,5,4,5,4,5,4,5,4,5,
5,4,4,4,4,4,6,5,6,4,5,4,5,4,5,4,4,5,6,5,4,4,4,
5,6,7,4,4,4,4,7,4,5,5,4,4,6,4,4,4,4,4,6,5,5,5,
4,3,6,4,4,5,4,4,5,4,5,4,4,5,4,4,3,4,4,8,7,4,4,
4,5,4,5,4,5,3,4,4,6,4,5,4,4,4,5,4,5,4,4,4,4,4,
4,5,4,5,4,6,6,4,5,4,4,5,6,5,6,4,4,4,5,4,4,6,7,
5,5,4,5,4,4,4,7,4,5,4,4,5,4,4,5,5,4,4,4,4,4,4,
5,6,4,4,7,5,5,4,7,5,4,4,4,6,4,4,5,6,4,5,5,5,4,
6,4,4,6,4,4,6,4,4,5,5,5,4,4,5,4,6,6,6,5,4,5,4,
5,6,5,5,6,6,6,5,4,3,4,4,5,6,5,4,5,4,4,4,4,4,5,
5,4,5,4,4,5,4,4,6,6,5,5,4,4,4,4,5,5,4,4,4,4,4,
4,4,4,5,6,4,5,4,4,4,4,5,4,4,4,5,5,4,4,4,4,4,4,
4,4,7,3,4,4,5,5,5,5,5,4,5,5,4,4,4,4,6,7,4,6,5,
4,4,6,4,5,5,4,6,5,5,5,5,5,4,6,5,4,5,4,5,5,6,6,
6,5,5,6,4,5,6,5,5,6,5,5,5,5,4,5,4,4,5,4,7,4,4,
4,4,4,7,5,6,5,4,4,5,6,7,4,6,5,4,4,6,5,5,5,4,4,
5,4,4,5,4,5,4,5,4,5,4,4,4,4,4,5,5,4,7,4,4,5,6,
4,5,6,4,4,3,5,4,4,4,4,4,4,7,5,5,6,3,4,4,5,6,4,
6,4,4,4,6,4,6,4,4,4,4,4,6,5,4,5,5,5,4,6,4,5,4,
4,4,4,4,3,4,5,6,4,3,6,4,4,6,4,5,5,4,4,5,4,4,5,
4,6,4,6,6,4,6,4,4,4,4,5,4,4,5,4,4,4,4,5,6,6,4,
6,4,4,4,4,4,3,6,4,5,4,6,4,5,2,7,4,4,5,5,4,6,4,
5,6,6,4,5,5,4,6,5,4,6,4,4,4,5,5,6,6,4,4,4,6,5,
5,4,5,6,4,4,6,7,4,4,4,4,4,5,6,4,4,4,4,5,4,4,4,
4,4,4,4,5,4,4,5,4,4,5,4,4,4,4,4,4,6,4,5,4,4,4,
4,5,5,4,4,4,5,4,5,5,5,4,4,6,5,4,4,5,4,4,4,5,4,
4,4,4,4,4,6,5,5,4,4,4,4,4,4,5,4,4,4,6,6,6,6,7,
5,4,6,3,5,5,4,5,4,3,5,4,6,6,4,4,5,4,4,4,5,4,5,
5,4,4,4,6,6,5,4,4,4,5,4,5,6,5,6,6,6,4,5,5,5,6,
4,5,6,5,4,4,5,4,5,4,4,6,4,4,4,4,5,4,4,4,4,4,4,
5,5,4,5,4,4,6,6,5,4,4,4,4,4,4,5,5,4,6,4,4,6,4,
7,7,5,5,5,5,5,5,4,4,4,4,5,4,5,6,4,6,5,5,4,6,6,
4,4,4,4,4,4,4,5,4,4,4,7,4,3,5,4,6,3,4,4,4,5,5,
4,6,6,6,4,5,5,6,5,6,6,5,6,6,6,4,5,5,6,4,5,5,5,
4,6,5,5,4,5,6,5,5,4,4,4,4,4,4,4,6,4,4,6,5,5,4,
4,4,4,6,6,5,4,4,6,5,4,4,4,5,6,6,4,6,4,4,5,6,4,
5,4,6,4,5,6,6,5,5,6,6,6,6,6,6,6,6,6,4,6,5,4,4,
6,8,6,7,6,6,6,6,6,6,6,6,6,6,6,7,6,6,6,6,7,6,6,
6,6,5,7,6,7,6,6,6,6,4,6,6,7,5,4,3,6,6,6,5,6,6,
6,6,6,6,6,6,6,6,6,6,6,5,6,6,6,6,6,7,4,5,5,4,5,
5,6,6,4,6,5,6,6,5,6,6,6,6,6,4,5,6,3,6,4,6,6,4,
5,5,6,6,4,5,6,4,5,5,4,6,4,5,5,4,6,6,6,3,6,6,5,
6,4,6,6,6,6,5,5,4,6,6,7,6,5,4,6,5,6,6,4,5,5,6,
4,6,6,6,6,5,6,5,7,5,4,5,6,5,6,5,6,6,6,6,5,4,6,
6,4,3,4,5,5,6,5,4,5,5,5,4,5,4,4,6,5,5,6,6,6,4,
5,5,6,5,3,6,6,6,3,5,6,6,4,4,4,4,5,4,4,4,6,5,4,
6,5,5,4,5,4,4,5,6,3,7,6,6,5,6,6,5,5,4,5,4,6,5,
4,5,6,4,5,4,3,4,4,6,5,6,6,4,6,4,5,4,5,4,6,6,5,
3,6,5,5,6,6,5,5,4,6,6,5,6,6,6,5,4,5,6,5,6,6,5,
6,5,5,4,5,4,6,6,5,6,4,5,5,6,5,5,6,5,5,5,5,5,6,
4,6,5,5,4,4,4,6,4,4,4,5,6,6,5,4,4,5,4,5,4,5,6,
6,5,4,5,4,5,5,6,4,4,4,4,4,4,4,6,5,5,5,5,4,4,4,
6,6,4,5,6,6,5,5,6,4,3,4,6,6,5,5,6,6,6,6,4,5,4,
4,4,6,4,5,3,3,5,4,6,6,6,6,6,4,6,4,4,4,6,4,5,4,
4,3,6,4,6,5,6,6,6,5,6,6,6,4,4,4,3,6,5,4,4,4,5,
4,4,4,5,4,4,5,4,5,3,4,4,3,4,6,4,4,5,4,7,4,4,5,
8,6,5,5,5,5,5,6,5,6,5,6,5,6,6,4,6,5,6,6,6,5,5,
5,5,5,5,6,6,5,5,5,5,6,5,6,5,6,6,5,6,5,5,5,5,5,
5,6,5,5,5,5,6,5,6,6,5,5,5,5,5,3,5,5,5,5,3,3,5,
5,5,5,5,4,6,5,5,5,4,5,5,5,4,5,5,4,5,5,4,5,4,4,
6,5,5,4,5,5,5,8,5,5,6,5,5,5,5,5,5,5,5,5,5,5,7,
5,5,6,6,6,5,5,6,6,5,5,5,5,5,6,6,6,5,6,5,6,5,5,
5,6,4,4,6,3,6,6,5,5,5,5,5,5,5,5,5,5,6,7,5,6,5,
5,4,4,5,5,5,5,6,5,5,6,5,5,5,4,5,5,5,4,4,5,4,4,
6,5,6,4,5,5,4,5,6,5,5,5,3,4,5,6,5,5,4,3,5,6,4,
4,5,4,5,5,3,6,4,5,4,5,5,4,5,3,3,5,5,5,5,5,6,5,
5,4,5,5,5,6,5,5,4,5,4,4,4,5,5,4,5,6,5,4,5,6,4,
6,5,5,3,5,5,4,5,5,6,7,6,5,4,5,4,5,5,4,5,5,4,4,
4,5,5,4,4,6,5,4,4,5,5,4,4,3,4,4,5,5,4,5,5,3,4,
4,4,6,6,6,5,4,4,4,3,6,5,5,4,6,5,5,4,4,4,6,5,6,
6,5,5,3,6,6,4,4,3,5,4,4,5,5,4,5,5,4,5,4,6,5,6,
4,4,4,4,5,6,5,4,4,4,6,4,4,4,6,5,4,4,5,4,5,3,4,
5,4,4,5,5,4,4,4,5,5,4,5,3,7,5,6,6,6,6,4,4,5,4,
4,4,5,4,4,5,4,5,4,6,4,5,4,5,5,6,4,3,3,5,5,6,5,
4,4,5,5,4,4,5,6,5,4,5,4,4,4,4,6,4,5,4,5,4,4,4,
4,4,5,4,5,4,6,4,3,3,6,5,5,6,6,6,5,6,4,4,4,5,6,
5,5,5,5,4,4,4,5,5,5,4,5,5,5,5,4,4,5,4,6,3,6,5,
5,3,5,6,5,6,6,5,6,5,4,4,5,4,4,4,4,5,4,5,4,4,5,
4,4,5,4,4,5,5,5,4,5,4,7,5,5,5,6,3,5,6,5,4,4,5,
4,5,5,4,6,4,4,4,4,4,5,4,4,5,5,4,4,5,6,6,4,6,5,
5,4,7,6,5,6,5,5,5,5,4,5,4,4,4,4,6,4,4,5,4,5,4,
5,5,5,4,4,4,5,5,5,4,4,5,4,6,5,5,5,6,5,5,5,5,5,
6,6,4,5,4,5,5,5,5,4,4,4,4,6,6,5,4,4,5,4,5,4,4,
4,4,5,6,4,4,4,3,5,5,4,5,6,4,4,5,2,4,5,4,5,5,6,
6,6,5,5,5,5,5,4,4,3,4,5,5,4,5,5,4,5,4,5,4,5,5,
5,4,4,5,4,4,4,4,4,4,4,5,3,3,4,4,4,5,7,5,6,3,6,
4,4,5,3,4,3,4,5,5,5,5,3,3,4,4,6,5,4,4,4,4,5,4,
4,5,4,4,4,6,5,5,4,4,4,4,4,4,4,4,4,4,4,4,5,5,4,
3,4,4,5,3,7,4,4,6,4,4,4,4,3,4,4,5,3,5,3,3,2,4,
4,4,4,3,4,5,4,4,4,4,5,4,3,3,5,5,5,4,6,4,5,5,4,
4,5,3,4,4,3,5,4,4,4,4,6,5,6,3,5,3,4,4,4,4,6,6,
5,6,5,4,4,4,4,5,6,5,4,4,5,4,5,3,4,4,4,3,4,5,5,
4,5,6,4,5,4,4,6,4,5,4,4,4,5,5,4,4,5,4,5,5,4,4,
4,4,4,4,5,4,2,5,4,5,6,5,6,5,7,4,4,3,5,5,4,4,4,
4,5,5,4,4,7,6,4,4,4,4,3,5,4,3,6,5,5,5,4,6,5,4,
4,4,5,5,4,5,4,5,4,6,4,4,4,4,5,4,5,5,4,4,5,4,5,
5,6,6,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,
4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,5,4,
4,4,4,4,4,4,6,4,4,4,4,4,4,4,4,4,4,4,4,4,5,6,6,
4,4,4,4,4,6,4,5,4,5,4,4,4,4,5,4,4,5,4,6,4,4,4,
4,4,4,4,4,5,5,5,4,5,5,4,4,4,4,5,5,4,5,4,4,4,4,
6,4,4,4,4,4,5,4,5,4,4,6,5,3,5,4,5,4,4,4,4,4,7,
5,5,5,5,4,3,4,5,4,4,5,3,4,4,5,4,3,4,4,4,5,4,4,
4,5,4,6,5,4,4,5,5,4,4,5,5,4,4,4,5,5,5,4,5,4,4,
5,5,3,4,4,4,6,4,5,6,5,4,4,4,5,6,5,5,4,5,4,5,4,
5,4,4,5,5,6,5,4,4,6,5,2,4,4,4,3,4,5,6,5,4,4,4,
4,4,4,4,4,4,4,4,5,5,5,5,4,5,4,5,6,4,5,4,5,5,5,
6,5,5,4,4,5,4,4,7,4,5,5,5,4,4,4,4,6,5,4,5,4,4,
4,3,5,4,4,5,4,6,4,5,4,5,4,5,5,5,5,5,6,5,4,5,4,
4,5,4,4,4,5,5,5,5,5,4,5,4,4,4,5,5,4,4,4,4,6,5,
6,4,4,4,4,6,6,6,6,5,4,6,5,4,4,4,5,4,6,5,4,6,4,
5,5,5,6,6,4,4,5,5,4,5,5,4,5,5,5,4,3,5,5,4,4,6,
7,5,5,5,5,5,5,5,5,5,5,6,5,5,5,5,5,5,5,6,6,4,5,
5,5,6,6,5,6,3,4,5,4,5,5,5,5,5,4,6,4,5,4,4,4,4,
4,4,4,4,4,4,5,4,4,4,4,4,4,5,5,4,4,5,4,4,4,5,5,
6,5,4,4,5,4,4,5,4,5,5,6,5,4,6,5,6,6,6,4,5,5,6,
5,5,6,5,5,5,5,4,4,5,5,5,5,4,4,5,4,4,4,5,5,4,7,
5,5,5,6,4,4,4,6,4,4,4,7,5,5,4,4,5,5,4,6,6,4,5,
5,5,5,4,5,4,7,6,5,5,5,5,5,4,5,5,5,4,5,4,4,5,5,
4,4,4,4,4,4,4,5,4,4,5,4,3,5,5,4,5,4,4,4,5,4,4,
6,5,5,4,5,5,5,4,4,4,4,5,5,4,4,6,4,4,5,6,6,4,4,
5,4,3,7,6,5,6,6,5,5,5,5,5,7,4,5,4,4,4,4,4,4,4,
4,5,4,4,4,4,5,4,5,5,5,4
],
  show_simple_stats(Magnitudes),
  show_histogram(Magnitudes),
  show_hpd_intervals(Magnitudes,[0.84,0.94,0.99]),
  [Methods,Explain] = recommend_abc_stats_explained(Magnitudes),
  println(methods=Methods),
  println(explain=Explain),
  reset_store,
  run_model(10_000,$model(Magnitudes),[show_probs_trunc,mean,
                           % show_percentiles,
                           show_hpd_intervals,hpd_intervals=[0.84,0.94,0.99],
                           % show_histogram,
                           min_accepted_samples=100,show_accepted_samples=true,
                           garbage_collect_when=10
                          ]),
  nl,
  println("Post:"),
  Post = get_store().get("post"),
  show_simple_stats(Post),
  show_histogram(Post),
  show_hpd_intervals(Post,[0.84,0.94,0.99]),  
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model(Magnitudes) =>
  K = uniform(0.1,10),
  Alpha = uniform(0.1,10),
  
  Y = pareto1_dist_n(K,Alpha,Magnitudes.len),
  observe_abc(Magnitudes,Y,1,[mean,stdev,q25,q75]),

  Post = pareto1_dist(K,Alpha),
  P = check(Post >= 6),

  if observed_ok then
    add("k",K),
    add("alpha",Alpha),
    add("post",Post),
    add("p",P)
  end.