/* Poisson distribution recover parameters in Picat. Generating random (poisson 20) and trying to recover data Cf my Gamble model gamble_poisson_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. /* Num accepted samples: 1000 Total samples: 3072 (0.326%) var : lambda Probabilities (truncated): 21.031627112752659: 0.0010000000000000 21.006523964255123: 0.0010000000000000 20.955794016503241: 0.0010000000000000 20.93465894118825: 0.0010000000000000 ......... 19.155100890600409: 0.0010000000000000 19.125061818690234: 0.0010000000000000 19.103388250140949: 0.0010000000000000 18.855326020539088: 0.0010000000000000 mean = 19.987 [min = 18.8553,mean = 19.987,median = 19.9765,max = 21.0316,variance = 0.116795,stdev = 0.341754] HPD intervals: HPD interval (0.84): 19.55152035241400..20.52393691130705 HPD interval (0.9): 19.45587485148741..20.56306760718263 HPD interval (0.95): 19.30769055888996..20.59816243213648 HPD interval (0.99): 19.26398326592652..20.86756951781473 Percentiles: (0.001 19.103140187911347) (0.01 19.292660560182274) (0.025 19.350627349413163) (0.05 19.436599903276118) (0.1 19.553523108441382) (0.25 19.758066845936849) (0.5 19.976487967224891) (0.75 20.220760345757789) (0.84 20.324303740150274) (0.9 20.44583232679528) (0.95 20.552762406527673) (0.975 20.694878406508863) (0.99 20.812146965475769) (0.999 21.006549067403622) (0.9999 21.029119308217759) (0.99999 21.031376332299171) Histogram (total 1000) 18.855: 1 # (0.001 / 0.001) 18.910: 0 (0.000 / 0.001) 18.964: 0 (0.000 / 0.001) 19.019: 0 (0.000 / 0.001) 19.073: 0 (0.000 / 0.001) 19.127: 2 ## (0.002 / 0.003) 19.182: 2 ## (0.002 / 0.005) 19.236: 1 # (0.001 / 0.006) 19.291: 10 ######## (0.010 / 0.016) 19.345: 16 ############## (0.016 / 0.032) 19.399: 18 ############### (0.018 / 0.050) 19.454: 20 ################# (0.020 / 0.070) 19.508: 23 ################### (0.023 / 0.093) 19.563: 35 ############################## (0.035 / 0.128) 19.617: 31 ########################## (0.031 / 0.159) 19.671: 37 ############################### (0.037 / 0.196) 19.726: 48 ######################################### (0.048 / 0.244) 19.780: 66 ######################################################## (0.066 / 0.310) 19.835: 66 ######################################################## (0.066 / 0.376) 19.889: 55 ############################################## (0.055 / 0.431) 19.943: 60 ################################################### (0.060 / 0.491) 19.998: 71 ############################################################ (0.071 / 0.562) 20.052: 62 #################################################### (0.062 / 0.624) 20.107: 40 ################################## (0.040 / 0.664) 20.161: 53 ############################################# (0.053 / 0.717) 20.216: 56 ############################################### (0.056 / 0.773) 20.270: 47 ######################################## (0.047 / 0.820) 20.324: 34 ############################# (0.034 / 0.854) 20.379: 25 ##################### (0.025 / 0.879) 20.433: 29 ######################### (0.029 / 0.908) 20.488: 25 ##################### (0.025 / 0.933) 20.542: 23 ################### (0.023 / 0.956) 20.596: 8 ####### (0.008 / 0.964) 20.651: 8 ####### (0.008 / 0.972) 20.705: 9 ######## (0.009 / 0.981) 20.760: 5 #### (0.005 / 0.986) 20.814: 7 ###### (0.007 / 0.993) 20.868: 3 ### (0.003 / 0.996) 20.923: 1 # (0.001 / 0.997) 20.977: 3 ### (0.003 / 1.000) var : post Probabilities (truncated): 20: 0.0960000000000000 21: 0.0880000000000000 19: 0.0850000000000000 17: 0.0810000000000000 ......... 36: 0.0020000000000000 10: 0.0020000000000000 34: 0.0010000000000000 5: 0.0010000000000000 mean = 20.125 [min = 5,mean = 20.125,median = 20.0,max = 37,variance = 19.8134,stdev = 4.45122] HPD intervals: HPD interval (0.84): 13.00000000000000..25.00000000000000 HPD interval (0.9): 13.00000000000000..26.00000000000000 HPD interval (0.95): 11.00000000000000..27.00000000000000 HPD interval (0.99): 10.00000000000000..32.00000000000000 Percentiles: (0.001 9.995000000000001) (0.01 11.99) (0.025 13) (0.05 13) (0.1 15) (0.25 17) (0.5 20) (0.75 23) (0.84 25) (0.9 26) (0.95 27) (0.975 30) (0.99 32) (0.999 37) (0.9999 37) (0.99999 37) Histogram (total 1000) 5: 1 # (0.001 / 0.001) 10: 2 # (0.002 / 0.003) 11: 7 #### (0.007 / 0.010) 12: 11 ####### (0.011 / 0.021) 13: 30 ################### (0.030 / 0.051) 14: 44 ############################ (0.044 / 0.095) 15: 53 ################################# (0.053 / 0.148) 16: 76 ################################################ (0.076 / 0.224) 17: 81 ################################################### (0.081 / 0.305) 18: 72 ############################################# (0.072 / 0.377) 19: 85 ##################################################### (0.085 / 0.462) 20: 96 ############################################################ (0.096 / 0.558) 21: 88 ####################################################### (0.088 / 0.646) 22: 75 ############################################### (0.075 / 0.721) 23: 51 ################################ (0.051 / 0.772) 24: 65 ######################################### (0.065 / 0.837) 25: 47 ############################# (0.047 / 0.884) 26: 40 ######################### (0.040 / 0.924) 27: 29 ################## (0.029 / 0.953) 28: 6 #### (0.006 / 0.959) 29: 15 ######### (0.015 / 0.974) 30: 6 #### (0.006 / 0.980) 31: 9 ###### (0.009 / 0.989) 32: 3 ## (0.003 / 0.992) 33: 3 ## (0.003 / 0.995) 34: 1 # (0.001 / 0.996) 36: 2 # (0.002 / 0.998) 37: 2 # (0.002 / 1.000) var : p Probabilities: 0: 0.5580000000000001 1: 0.4420000000000000 mean = 0.442 [min = 0,mean = 0.442,median = 0.0,max = 1,variance = 0.246636,stdev = 0.496625] HPD intervals: HPD interval (0.84): 0.00000000000000..1.00000000000000 HPD interval (0.9): 0.00000000000000..1.00000000000000 HPD interval (0.95): 0.00000000000000..1.00000000000000 HPD interval (0.99): 0.00000000000000..1.00000000000000 Percentiles: (0.001 0) (0.01 0) (0.025 0) (0.05 0) (0.1 0) (0.25 0) (0.5 0) (0.75 1) (0.84 1) (0.9 1) (0.95 1) (0.975 1) (0.99 1) (0.999 1) (0.9999 1) (0.99999 1) Histogram (total 1000) 0: 558 ############################################################ (0.558 / 0.558) 1: 442 ################################################ (0.442 / 1.000) * Mean: Picat> X = poisson_dist_mean(20) X = 20 * Probablity of > 20: Picat> X = 1-poisson_dist_cdf(20,20) X = 0.440907415768675 */ go ?=> NumSamples = 250, Data = poisson_dist_n(20,NumSamples), println(data=Data), show_simple_stats(Data), reset_store, run_model(10_000,$model(Data),[show_probs_trunc,mean, show_simple_stats, show_hpd_intervals,hpd_intervals=[0.84,0.9,0.95,0.99], show_percentiles, show_histogram, min_accepted_samples=1000,show_accepted_samples=true ]), nl, % show_store_lengths,nl, % fail, nl. go => true. model(Data) => Lambda = normal_dist(Data.mean,1), X = poisson_dist_n(Lambda,Data.len), % observe_abc(Data,X), observe_abc(Data,X,1/10), Post = poisson_dist(Lambda), P = check1(Post > 20), if observed_ok then add("lambda",Lambda), add("post",Post), add("p",P), end.