/* 

  Sum Pareto in Picat.

  From PSI test/approximate/sumpareto.psi
  """
  def main() (
    s:=0;
    for i in (0..2) (
        x:=pareto(1.161,1); //80-20 rule: https://en.wikipedia.org/wiki/Pareto_distribution#Relation_to_the_.22Pareto_principle.22
		s = s + x;
    ) ;; mathematica couldn't solve this integral even after 5 days. 
    return s;
  )
  ;; E(s) = 2322/161
  ;; (14.4223602484472)
  """  

  Note:
  Picat> X = pareto1_dist_mean(1, 1.161)*2
  X = 14.422360248447202


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

main => go.

/*
  Here are some runs:

  var : s
  Probabilities (truncated):
  19925.267972070254473: 0.0000500000000000
  11595.132132638280382: 0.0000500000000000
  8241.184898646772126: 0.0000500000000000
  5252.708786347981913: 0.0000500000000000
  .........
  2.017744963785102: 0.0000500000000000
  2.017334971021555: 0.0000500000000000
  2.016290694521585: 0.0000500000000000
  2.01069967918622: 0.0000500000000000
  mean = 13.4563


  var : s
  Probabilities (truncated):
  6002.456595301281595: 0.0000500000000000
  2842.1291807196626: 0.0000500000000000
  2722.665508356536066: 0.0000500000000000
  2425.205637614298666: 0.0000500000000000
  .........
  2.01619013297561: 0.0000500000000000
  2.014021043210777: 0.0000500000000000
  2.012151656349488: 0.0000500000000000
  2.007954712923815: 0.0000500000000000
  mean = 11.1607

  var : s
  Probabilities (truncated):
  10531.169376848152751: 0.0000500000000000
  9250.108613417267406: 0.0000500000000000
  5003.043126168410708: 0.0000500000000000
  4894.318180659392965: 0.0000500000000000
  .........
  2.018947415665721: 0.0000500000000000
  2.018347271805737: 0.0000500000000000
  2.015019115683693: 0.0000500000000000
  2.005032482597694: 0.0000500000000000
  mean = 12.373


*/
go ?=>
  reset_store,
  run_model(20_000,$model,[show_probs_trunc,mean]),
  nl,
  % show_store_lengths,
  % fail,
  nl.
go => true.
  
  
model() =>

  D = pareto1_dist_n(1,1.161, 2),
  S = D.sum,
  % S = 0,
  % foreach(_ in 1..2)
  %   S := S + pareto1_dist(1,1.161)
  % end,

  % add("f",F),
  add("s",S).