/* Geometric distribution in Picat. See ppl_distributions.pi and ppl_distributions_test.pi for more on this. Theoretical: * Mean: Picat> X=geometric_dist_mean(0.1) X = 9.0 * Probability of larger than 15 Picat> X=1-geometric_dist_cdf(0.1,15) X = 0.185302018885184 * Exactly 15 Picat> X=geometric_dist_pdf(0.1,15) X = 0.020589113209465 Picat> pdf_all($geometric_dist(1/10)).printf_list 0 0.1 1 0.09 2 0.081 3 0.0729 4 0.06561 5 0.059049 6 0.0531441 7 0.04782969 8 0.043046721 9 0.0387420489 10 0.03486784401 11 0.031381059609 12 0.0282429536481 13 0.02541865828329 14 0.022876792454961 15 0.020589113209465 16 0.018530201888518 17 0.016677181699667 18 0.0150094635297 19 0.01350851717673 20 0.012157665459057 21 0.010941898913151 22 0.009847709021836 23 0.008862938119653 24 0.007976644307687 25 0.007178979876919 26 0.006461081889227 27 0.005814973700304 28 0.005233476330274 29 0.004710128697246 30 0.004239115827522 31 0.003815204244769 32 0.003433683820293 33 0.003090315438263 34 0.002781283894437 35 0.002503155504993 36 0.002252839954494 37 0.002027555959045 38 0.00182480036314 39 0.001642320326826 40 0.001478088294143 41 0.001330279464729 42 0.001197251518256 43 0.001077526366431 Cf my Gamble model gamble_geometric_dist.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. /* var : d Probabilities: 0: 0.1065000000000000 1: 0.0872000000000000 2: 0.0850000000000000 3: 0.0711000000000000 4: 0.0676000000000000 5: 0.0574000000000000 6: 0.0489000000000000 7: 0.0468000000000000 8: 0.0416000000000000 9: 0.0377000000000000 10: 0.0358000000000000 11: 0.0346000000000000 12: 0.0273000000000000 14: 0.0238000000000000 13: 0.0229000000000000 16: 0.0203000000000000 15: 0.0198000000000000 17: 0.0173000000000000 18: 0.0164000000000000 20: 0.0141000000000000 19: 0.0132000000000000 21: 0.0123000000000000 22: 0.0097000000000000 23: 0.0086000000000000 24: 0.0071000000000000 26: 0.0065000000000000 25: 0.0060000000000000 29: 0.0057000000000000 28: 0.0051000000000000 27: 0.0049000000000000 31: 0.0037000000000000 32: 0.0033000000000000 30: 0.0033000000000000 33: 0.0032000000000000 35: 0.0028000000000000 34: 0.0022000000000000 36: 0.0020000000000000 37: 0.0019000000000000 40: 0.0016000000000000 43: 0.0014000000000000 42: 0.0014000000000000 41: 0.0013000000000000 38: 0.0013000000000000 39: 0.0010000000000000 50: 0.0009000000000000 44: 0.0007000000000000 54: 0.0005000000000000 49: 0.0005000000000000 45: 0.0005000000000000 62: 0.0004000000000000 58: 0.0004000000000000 52: 0.0004000000000000 48: 0.0004000000000000 82: 0.0003000000000000 59: 0.0003000000000000 57: 0.0003000000000000 53: 0.0003000000000000 47: 0.0003000000000000 46: 0.0003000000000000 66: 0.0002000000000000 65: 0.0002000000000000 61: 0.0002000000000000 60: 0.0002000000000000 56: 0.0002000000000000 96: 0.0001000000000000 76: 0.0001000000000000 74: 0.0001000000000000 72: 0.0001000000000000 71: 0.0001000000000000 67: 0.0001000000000000 63: 0.0001000000000000 55: 0.0001000000000000 51: 0.0001000000000000 mean = 8.8911 var : p Probabilities: false: 0.8139999999999999 true: 0.1860000000000000 mean = [false = 0.814,true = 0.186] var : p2 Probabilities: false: 0.9802000000000000 true: 0.0198000000000000 mean = [false = 0.9802,true = 0.0198] */ go ?=> reset_store, run_model(10_000,$model,[show_probs,mean]), nl, % show_store_lengths,nl, % fail, nl. go => true. model() => D = geometric_dist(1/10), P = check(D > 15), P2 = check(D == 15), add("d",D), add("p",P), add("p2",P2). /* Recover parameter from geometric_dist(1/10). ... Num accepted samples: 9997 Total samples: 218386 (0.046%) Num accepted samples: 9998 Total samples: 218397 (0.046%) Num accepted samples: 9999 Total samples: 218417 (0.046%) Num accepted samples: 10000 Total samples: 218437 (0.046%) var : p Probabilities (truncated): 0.186821053822907: 0.0001000000000000 0.185954484709517: 0.0001000000000000 0.18226094878384: 0.0001000000000000 0.182040868877452: 0.0001000000000000 ......... 0.073873618186393: 0.0001000000000000 0.073788928833692: 0.0001000000000000 0.073585193638497: 0.0001000000000000 0.07163865588216: 0.0001000000000000 mean = 0.116117 */ go2 ?=> % Data = geometric_dist_n(1/10,100), Data = [14,20,9,21,47,3,5,3,8,2,8,6,0,4,5,7,15,1,17,1,2,9,32,3,11,1,2,1,15,3,2,8,2,20,6,2,2,23,6,14,1,0,3,2,6,12,14,3,19,9,5,0,2,4,4,35,6,9,1,3,31,4,21,1,3,4,4,6,3,19,3,5,32,7,10,5,3,5,12,1,0,1,2,3,7,9,3,0,2,5,4,2,18,3,5,1,9,14,9,26], println(data=Data), println([mean=Data.mean,stdev=Data.stdev]), reset_store, run_model(10_000,$model2(Data),[show_probs_trunc,mean, min_accepted_samples=10_000,show_accepted_samples=true ]), nl, % show_store_lengths,nl, % fail, nl. go2 => true. model2(Data) => Len = Data.len, Mean = Data.mean, Stdev = Data.stdev, P = uniform(0,1), Y = geometric_dist_n(P,Len), % YMean = Y.mean, % YStdev = Y.stdev, % % Note the reduced Stdev % observe(abs(Mean-YMean)