/* 

  Adding to 7 in Picat.

  From  https://www.youtube.com/watch?v=6Lqt07enBGs&list=PLPhJUsEbNh6T7TZ-LPRZgepphZQIB0A02&index=6&t=1627s
  "Probabilistic Programming and Bayesian Nonparametrics -- Frank Wood (Part 1)"
  @1:21:01
  """  
   [assume a (- (poisson 100) 100)]
   [assume b (- (poisson 100) 100)]
   [observe (normal (+ a b) 0.00001) 7]
   [predict (list a b)]
  """

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

/*

var : a
Probabilities (truncated):
6: 0.0658064516129032
5: 0.0606451612903226
3: 0.0541935483870968
9: 0.0516129032258065
.........
-17: 0.0012903225806452
-18: 0.0012903225806452
-19: 0.0012903225806452
-22: 0.0012903225806452
mean = 3.32645
Histogram:
-22.000:     1 #                                                            (0.001 / 0.001)
-20.850:     0                                                              (0.000 / 0.001)
-19.700:     0                                                              (0.000 / 0.001)
-18.550:     2 #                                                            (0.003 / 0.004)
-17.400:     1 #                                                            (0.001 / 0.005)
-16.250:     3 ##                                                           (0.004 / 0.009)
-15.100:     2 #                                                            (0.003 / 0.012)
-13.950:     2 #                                                            (0.003 / 0.014)
-12.800:     4 ###                                                          (0.005 / 0.019)
-11.650:     4 ###                                                          (0.005 / 0.025)
-10.500:    18 ############                                                 (0.023 / 0.048)
 -9.350:    12 ########                                                     (0.015 / 0.063)
 -8.200:    11 ########                                                     (0.014 / 0.077)
 -7.050:    16 ###########                                                  (0.021 / 0.098)
 -5.900:    20 ##############                                               (0.026 / 0.124)
 -4.750:    11 ########                                                     (0.014 / 0.138)
 -3.600:    24 #################                                            (0.031 / 0.169)
 -2.450:    72 ##################################################           (0.093 / 0.262)
 -1.300:    34 #######################                                      (0.044 / 0.306)
 -0.150:    33 #######################                                      (0.043 / 0.348)
  1.000:    32 ######################                                       (0.041 / 0.390)
  2.150:    38 ##########################                                   (0.049 / 0.439)
  3.300:    42 #############################                                (0.054 / 0.493)
  4.450:    87 ############################################################ (0.112 / 0.605)
  5.600:    51 ###################################                          (0.066 / 0.671)
  6.750:    30 #####################                                        (0.039 / 0.710)
  7.900:    26 ##################                                           (0.034 / 0.743)
  9.050:    40 ############################                                 (0.052 / 0.795)
 10.200:    32 ######################                                       (0.041 / 0.836)
 11.350:    13 #########                                                    (0.017 / 0.853)
 12.500:    48 #################################                            (0.062 / 0.915)
 13.650:    22 ###############                                              (0.028 / 0.943)
 14.800:    12 ########                                                     (0.015 / 0.959)
 15.950:     6 ####                                                         (0.008 / 0.966)
 17.100:     9 ######                                                       (0.012 / 0.978)
 18.250:     4 ###                                                          (0.005 / 0.983)
 19.400:     4 ###                                                          (0.005 / 0.988)
 20.550:     8 ######                                                       (0.010 / 0.999)
 21.700:     0                                                              (0.000 / 0.999)
 22.850:     1 #                                                            (0.001 / 1.000)

var : b
Probabilities (truncated):
6: 0.0658064516129032
5: 0.0606451612903226
3: 0.0541935483870968
9: 0.0516129032258065
.........
-17: 0.0012903225806452
-18: 0.0012903225806452
-19: 0.0012903225806452
-22: 0.0012903225806452
mean = 3.32645
Histogram:
-22.000:     1 #                                                            (0.001 / 0.001)
-20.850:     0                                                              (0.000 / 0.001)
-19.700:     0                                                              (0.000 / 0.001)
-18.550:     2 #                                                            (0.003 / 0.004)
-17.400:     1 #                                                            (0.001 / 0.005)
-16.250:     3 ##                                                           (0.004 / 0.009)
-15.100:     2 #                                                            (0.003 / 0.012)
-13.950:     2 #                                                            (0.003 / 0.014)
-12.800:     4 ###                                                          (0.005 / 0.019)
-11.650:     4 ###                                                          (0.005 / 0.025)
-10.500:    18 ############                                                 (0.023 / 0.048)
 -9.350:    12 ########                                                     (0.015 / 0.063)
 -8.200:    11 ########                                                     (0.014 / 0.077)
 -7.050:    16 ###########                                                  (0.021 / 0.098)
 -5.900:    20 ##############                                               (0.026 / 0.124)
 -4.750:    11 ########                                                     (0.014 / 0.138)
 -3.600:    24 #################                                            (0.031 / 0.169)
 -2.450:    72 ##################################################           (0.093 / 0.262)
 -1.300:    34 #######################                                      (0.044 / 0.306)
 -0.150:    33 #######################                                      (0.043 / 0.348)
  1.000:    32 ######################                                       (0.041 / 0.390)
  2.150:    38 ##########################                                   (0.049 / 0.439)
  3.300:    42 #############################                                (0.054 / 0.493)
  4.450:    87 ############################################################ (0.112 / 0.605)
  5.600:    51 ###################################                          (0.066 / 0.671)
  6.750:    30 #####################                                        (0.039 / 0.710)
  7.900:    26 ##################                                           (0.034 / 0.743)
  9.050:    40 ############################                                 (0.052 / 0.795)
 10.200:    32 ######################                                       (0.041 / 0.836)
 11.350:    13 #########                                                    (0.017 / 0.853)
 12.500:    48 #################################                            (0.062 / 0.915)
 13.650:    22 ###############                                              (0.028 / 0.943)
 14.800:    12 ########                                                     (0.015 / 0.959)
 15.950:     6 ####                                                         (0.008 / 0.966)
 17.100:     9 ######                                                       (0.012 / 0.978)
 18.250:     4 ###                                                          (0.005 / 0.983)
 19.400:     4 ###                                                          (0.005 / 0.988)
 20.550:     8 ######                                                       (0.010 / 0.999)
 21.700:     0                                                              (0.000 / 0.999)
 22.850:     1 #                                                            (0.001 / 1.000)

var : a+b
Probabilities:
7: 1.0000000000000000
mean = 7.0
Histogram:
      7:   775 ############################################################ (1.000 / 1.000)

*/
go ?=>
  reset_store,
  run_model(30_000,$model,[show_probs_trunc,mean,show_histogram]),
  nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  N = 100,

  A = poisson_dist(N) - N,
  B = poisson_dist(N) - N,  
  
  observe( abs(  normal_dist(A+B,0.0001)-7) <= 0.01),
  
  if observed_ok then
    add("a",A),
    add("b",A),
    add("a+b",A+B)
  end.