/* 

  Zip dist (recover parameters) in Picat.

  Recover parameters for zipf 15 0.3

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

/*

  Not very bad, but not good either.

data = [12,5,11,12,14,2,4,11,3,7,6,9,5,7,15,14,9,10,2,8,1,3,2,12,2,5,2,1,15,3,7,12,8,4,9,7,6,15,4,11,7,11,5,13,3,4,12,14,1,14,7,1,2,9,13,4,1,1,6,1,3,15,14,13,3,7,5,11,7,9,7,1,6,14,14,10,3,11,9,4,10,2,6,13,12,4,3,13,4,10,15,8,9,13,6,14,5,12,10,14]
[min = 1,mean = 7.73,median = 7.0,max = 15,variance = 19.6371,stdev = 4.43138]
Histogram (total 100)
      1:     8 #####################################################        (0.080 / 0.080)
      2:     7 ###############################################              (0.070 / 0.150)
      3:     8 #####################################################        (0.080 / 0.230)
      4:     8 #####################################################        (0.080 / 0.310)
      5:     6 ########################################                     (0.060 / 0.370)
      6:     6 ########################################                     (0.060 / 0.430)
      7:     9 ############################################################ (0.090 / 0.520)
      8:     3 ####################                                         (0.030 / 0.550)
      9:     7 ###############################################              (0.070 / 0.620)
     10:     5 #################################                            (0.050 / 0.670)
     11:     6 ########################################                     (0.060 / 0.730)
     12:     7 ###############################################              (0.070 / 0.800)
     13:     6 ########################################                     (0.060 / 0.860)
     14:     9 ############################################################ (0.090 / 0.950)
     15:     5 #################################                            (0.050 / 1.000)
theoretical_mean = 7.01817
methods = [mean,stdev,q25,q75]
explain = Data shows mild skewness or moderate tails. Adding quantiles for shape.
min_accepted_samples = 100
var : n
Probabilities:
15: 0.4000000000000000
16: 0.3000000000000000
14: 0.1800000000000000
17: 0.1200000000000000
mean = 15.36

var : s
Probabilities (truncated):
0.51999703073874: 0.0100000000000000
0.483201695831121: 0.0100000000000000
0.45955737142803: 0.0100000000000000
0.458897850689896: 0.0100000000000000
.........
0.021977135921818: 0.0100000000000000
0.015199109918996: 0.0100000000000000
0.009497082796645: 0.0100000000000000
0.001021316275476: 0.0100000000000000
mean = 0.188839

var : post
Probabilities (truncated):
1: 0.1200000000000000
7: 0.1000000000000000
6: 0.1000000000000000
3: 0.1000000000000000
.........
12: 0.0400000000000000
4: 0.0400000000000000
15: 0.0300000000000000
10: 0.0300000000000000
mean = 7.07


Post:
[min = 1,mean = 7.07,median = 7.0,max = 15,variance = 18.3251,stdev = 4.28078]
Histogram (total 100)
      1:    12 ############################################################ (0.120 / 0.120)
      2:     6 ##############################                               (0.060 / 0.180)
      3:    10 ##################################################           (0.100 / 0.280)
      4:     4 ####################                                         (0.040 / 0.320)
      5:     6 ##############################                               (0.060 / 0.380)
      6:    10 ##################################################           (0.100 / 0.480)
      7:    10 ##################################################           (0.100 / 0.580)
      8:     6 ##############################                               (0.060 / 0.640)
      9:     7 ###################################                          (0.070 / 0.710)
     10:     3 ###############                                              (0.030 / 0.740)
     11:     6 ##############################                               (0.060 / 0.800)
     12:     4 ####################                                         (0.040 / 0.840)
     13:     5 #########################                                    (0.050 / 0.890)
     14:     8 ########################################                     (0.080 / 0.970)
     15:     3 ###############                                              (0.030 / 1.000)


*/
go ?=>
  Data = zipf_dist_n(15,0.3,100),
  println(data=Data),
  show_simple_stats(Data),
  Data.show_histogram(),
  println(theoretical_mean=zipf_dist_mean(15,0.3)),
  [Methods,Explain] = recommend_abc_stats_explained(Data),
  println(methods=Methods),
  println(explain=Explain),
  reset_store,
  run_model(10_000,$model(Data),[show_probs_trunc,mean,
                           % show_percentiles,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % show_histogram,
                           min_accepted_samples=100,show_accepted_samples=false % true
                          ]),
  nl,
  println("Post:"),
  Post = get_store().get("post"),
  show_simple_stats(Post),
  Post.show_histogram(),
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model(Data) =>
  N = random_integer1(2*Data.max),
  S = uniform(0,10),

  Y = zipf_dist_n(N,S,Data.len),
  observe_abc(Data,Y,1/10),
  % observe_abc(Data,Y,1/2,[mean,stdev,q25,q75]), % worse
  
  Post = zipf_dist(N,S),

  if observed_ok then
    add("n",N),
    add("s",S),
    add("post",Post)    
  end.