/* Exponential dist in Picat. From Handbook on probability distributions page 56ff See ppl_distributions.pi and ppl_distributions_test.pi for more on this, e.g. PDF, CDF, quantile, mean, and variance. From Mathematica's ExponentialDistribution """ A battery has a lifespan that is exponentially distributed with rate parameter 1/3000 per hour. Find the probability that a random battery has a lifespan of less than 2500 hours: Probability[t < 2500.0, t e ExponentialDistribution[1/3000]] -> 0.565402 """ Picat> X=exponential_dist_pdf(1/3000, 2500) X = 0.000144866069502 CDF and Quantile: Picat> X=exponential_dist_cdf(1/3000, 2500), exponential_dist_quantile(1/3000,X)=Y X = 0.565401791492922 Y = 2499.999999999999545 Cf my Gamble model gamble_.exponential_distrkt 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 : d1 Probabilities (truncated): 33312.456075991023681: 0.0001000000000000 30286.31540420496458: 0.0001000000000000 23856.428990803011402: 0.0001000000000000 23183.477646950257622: 0.0001000000000000 ......... 2.117339614015015: 0.0001000000000000 2.038570620710133: 0.0001000000000000 1.16659795925845: 0.0001000000000000 0.742854314389566: 0.0001000000000000 mean = 3011.95 Percentiles: (0.001 4.8791945446479952) (0.01 30.296083778756572) (0.025 83.311120510646774) (0.05 160.71056527862254) (0.1 306.93017373055483) (0.25 861.54486008674121) (0.5 2090.8780846834579) (0.75 4139.547615119578) (0.84 5537.0085743396276) (0.9 6861.4856802174063) (0.95 9030.5188565104054) (0.975 11368.900541547591) (0.99 14274.42209813129) (0.999 21944.540594100934) (0.9999 30286.618018270005) (0.99999 33009.872270222775) HPD intervals: HPD interval (0.84): 1.16659795925845..5537.15376731525248 Histogram (total 10000) 0.743: 1300 ##################################### (0.130 / 0.130) 833.536: 2109 ############################################################ (0.211 / 0.341) 1666.329: 1584 ############################################# (0.158 / 0.499) 2499.121: 1220 ################################### (0.122 / 0.621) 3331.914: 916 ########################## (0.092 / 0.713) 4164.707: 720 #################### (0.072 / 0.785) 4997.500: 483 ############## (0.048 / 0.833) 5830.293: 434 ############ (0.043 / 0.877) 6663.085: 316 ######### (0.032 / 0.908) 7495.878: 212 ###### (0.021 / 0.929) 8328.671: 168 ##### (0.017 / 0.946) 9161.464: 105 ### (0.011 / 0.957) 9994.257: 103 ### (0.010 / 0.967) 10827.050: 73 ## (0.007 / 0.974) 11659.842: 53 ## (0.005 / 0.980) 12492.635: 56 ## (0.006 / 0.985) 13325.428: 35 # (0.004 / 0.989) 14158.221: 23 # (0.002 / 0.991) 14991.014: 18 # (0.002 / 0.993) 15823.807: 15 (0.002 / 0.994) 16656.599: 16 (0.002 / 0.996) 17489.392: 10 (0.001 / 0.997) 18322.185: 7 (0.001 / 0.998) 19154.978: 4 (0.000 / 0.998) 19987.771: 5 (0.001 / 0.998) 20820.564: 1 (0.000 / 0.999) 21653.356: 5 (0.001 / 0.999) 22486.149: 4 (0.000 / 0.999) 23318.942: 2 (0.000 / 1.000) 24151.735: 1 (0.000 / 1.000) 24984.528: 0 (0.000 / 1.000) 25817.321: 0 (0.000 / 1.000) 26650.113: 0 (0.000 / 1.000) 27482.906: 0 (0.000 / 1.000) 28315.699: 0 (0.000 / 1.000) 29148.492: 0 (0.000 / 1.000) 29981.285: 1 (0.000 / 1.000) 30814.078: 0 (0.000 / 1.000) 31646.870: 0 (0.000 / 1.000) 32479.663: 1 (0.000 / 1.000) var : p1 Probabilities: true: 0.5621000000000000 false: 0.4379000000000000 mean = [true = 0.5621,false = 0.4379] Percentiles: (0.001 false) (0.01 false) (0.025 false) (0.05 false) (0.1 false) (0.25 false) (0.5 true) (0.75 true) (0.84 true) (0.9 true) (0.95 true) (0.975 true) (0.99 true) (0.999 true) (0.9999 true) (0.99999 true) HPD intervals: show_hpd_intervals: data is not numeric Histogram (total 10000) false : 4379 ############################################### (0.438 / 0.438) true : 5621 ############################################################ (0.562 / 1.000) */ go ?=> reset_store, run_model(10_000,$model,[show_probs_trunc,mean,show_percentiles, show_hpd_intervals,hpd_intervals=[0.84], show_histogram]), nl, % show_store_lengths,nl, % fail, nl. go => true. model() => Lambda = 1/3000, % Note: Gamble has this as 3000 D1 = exponential_dist(Lambda), P1 = check(D1 < 2500), add("d1",D1), add("p1",P1). /* Recovering data: LambdaInv = 3000, exponential_dist(1/LambdaInv) * uniform(1,10000): Num accepted samples: 10000 Total samples: 394770 (0.025%) var : 1/lambda Probabilities (truncated): 7229.232217458650666: 0.0001000000000000 7126.098850006749672: 0.0001000000000000 6740.895767715245711: 0.0001000000000000 6636.740832252307882: 0.0001000000000000 ......... 1684.313267659076018: 0.0001000000000000 1679.644814248031253: 0.0001000000000000 1625.016455312732887: 0.0001000000000000 1503.862927443284207: 0.0001000000000000 mean = 3118.52 HPD intervals: HPD interval (0.84): 2240.96008111534638..3870.25007945031393 * random_integer1(10000) Num accepted samples: 10000 Total samples: 399793 (0.025%) var : 1/lambda Probabilities (truncated): 2833: 0.0015000000000000 2735: 0.0015000000000000 3126: 0.0014000000000000 3115: 0.0014000000000000 ......... 1667: 0.0001000000000000 1644: 0.0001000000000000 1639: 0.0001000000000000 1594: 0.0001000000000000 mean = 3113.16 HPD intervals: HPD interval (0.84): 2239.00000000000000..3899.00000000000000 */ go2 ?=> LambdaInv = 3000, Data = exponential_dist_n(1/LambdaInv,30), println([lambda=Lambda,data=Data]), reset_store, run_model(10_000,$model2(Data),[show_probs_trunc,mean, show_hpd_intervals,hpd_intervals=[0.84], min_accepted_samples=10000,show_accepted_samples=true ]), nl, % show_store_lengths,nl, % fail, nl. go2 => true. model2(Data) => % LambdaInv = uniform(1,10000), LambdaInv = random_integer1(10000), X = exponential_dist_n(1/LambdaInv,Data.len), observe_abc(Data,X,1/10), if observed_ok then add("1/lambda",LambdaInv) end.