/* Mean of uniform in Picat. How much differs the mean from random selection of (discrete) uniform distribution? Here we take n samples from 1..n. For n=10, there is a small chance that the mean is >= 8, about 0.0023 The exact value (from Model 2) is 0.0028197500000000115. * Model 1: Using randomInteger (1..10) with Rejection * Model 2: Using multinomial * Model 3: Parameter estimation from the results in Model 2 (100 samples) Cf my Gamble model gamble_mean_of_uniform.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. /* Using random_integer1 var : m Probabilities (truncated): 5.5: 0.0437000000000000 5.6: 0.0425000000000000 5.7: 0.0423000000000000 5.4: 0.0413000000000000 ......... 8.6: 0.0001000000000000 8.5: 0.0001000000000000 2.8: 0.0001000000000000 2.5: 0.0001000000000000 mean = 5.50387 Percentiles: (0.001 2.8999999999999999) (0.01 3.5) (0.025 3.7000000000000002) (0.05 4) (0.1 4.2999999999999998) (0.25 4.9000000000000004) (0.5 5.5) (0.75 6.0999999999999996) (0.84 6.4000000000000004) (0.9 6.7000000000000002) (0.95 7) (0.975 7.2999999999999998) (0.99 7.5999999999999996) (0.999 8.2001000000000204) (0.9999 8.6000399999997175) (0.99999 8.960004000000481) var : s Probabilities (truncated): 55: 0.0437000000000000 56: 0.0425000000000000 57: 0.0423000000000000 54: 0.0413000000000000 ......... 86: 0.0001000000000000 85: 0.0001000000000000 28: 0.0001000000000000 25: 0.0001000000000000 mean = 55.0387 Percentiles: (0.001 29) (0.01 35) (0.025 37) (0.05 40) (0.1 43) (0.25 49) (0.5 55) (0.75 61) (0.84 64) (0.9 67) (0.95 70) (0.975 73) (0.99 76) (0.999 82.001000000000204) (0.9999 86.000399999997171) (0.99999 89.60004000000481) var : p Probabilities: false: 0.9976000000000000 true: 0.0024000000000000 mean = 0.0024 Percentiles: (0.001 false) (0.01 false) (0.025 false) (0.05 false) (0.1 false) (0.25 false) (0.5 false) (0.75 false) (0.84 false) (0.9 false) (0.95 false) (0.975 false) (0.99 false) (0.999 true) (0.9999 true) (0.99999 true) */ go ?=> reset_store, run_model(10_000,$model,[show_probs_trunc,mean, show_percentiles ]), nl, % show_store_lengths,nl, % fail, nl. go => true. model() => N = 10, Sample = random_integer1_n(N,N), M = Sample.mean, S = Sample.sum, P = check(M >= 8), add("m",M), add("s",S), add("p",P). /* Using multinomial ar : m Probabilities (truncated): 5.7: 0.0423000000000000 5.2: 0.0422000000000000 5.5: 0.0420000000000000 5.6: 0.0418000000000000 ......... 2.5: 0.0002000000000000 8.699999999999999: 0.0001000000000000 8.5: 0.0001000000000000 2.0: 0.0001000000000000 mean = 5.5001 Percentiles: (0.001 2.7999999999999998) (0.01 3.3999999999999999) (0.025 3.7000000000000002) (0.05 4) (0.1 4.2999999999999998) (0.25 4.9000000000000004) (0.5 5.5) (0.75 6.0999999999999996) (0.84 6.4000000000000004) (0.9 6.7000000000000002) (0.95 7) (0.975 7.2999999999999998) (0.99 7.5999999999999996) (0.999 8.1999999999999993) (0.9999 8.6000099999999282) (0.99999 8.6900010000001195) Histogram (total 10000) 2.000: 1 (0.000 / 0.000) 2.167: 0 (0.000 / 0.000) 2.335: 0 (0.000 / 0.000) 2.502: 2 (0.000 / 0.000) 2.670: 5 (0.001 / 0.001) 2.837: 10 # (0.001 / 0.002) 3.005: 13 # (0.001 / 0.003) 3.172: 28 ## (0.003 / 0.006) 3.340: 71 ##### (0.007 / 0.013) 3.507: 50 #### (0.005 / 0.018) 3.675: 119 ######### (0.012 / 0.030) 3.842: 166 ############ (0.017 / 0.046) 4.010: 114 ######## (0.011 / 0.058) 4.178: 296 ##################### (0.030 / 0.087) 4.345: 415 ############################## (0.042 / 0.129) 4.512: 233 ################# (0.023 / 0.152) 4.680: 565 ######################################## (0.057 / 0.209) 4.848: 690 ################################################# (0.069 / 0.278) 5.015: 375 ########################### (0.037 / 0.315) 5.182: 836 ############################################################ (0.084 / 0.399) 5.350: 795 ######################################################### (0.080 / 0.478) 5.518: 838 ############################################################ (0.084 / 0.562) 5.685: 423 ############################## (0.042 / 0.605) 5.852: 770 ####################################################### (0.077 / 0.681) 6.020: 732 #################################################### (0.073 / 0.755) 6.188: 334 ######################## (0.033 / 0.788) 6.355: 598 ########################################### (0.060 / 0.848) 6.522: 443 ################################ (0.044 / 0.892) 6.690: 210 ############### (0.021 / 0.913) 6.857: 262 ################### (0.026 / 0.939) 7.025: 242 ################# (0.024 / 0.964) 7.192: 78 ###### (0.008 / 0.971) 7.360: 135 ########## (0.013 / 0.985) 7.527: 69 ##### (0.007 / 0.992) 7.695: 26 ## (0.003 / 0.994) 7.862: 26 ## (0.003 / 0.997) 8.030: 18 # (0.002 / 0.999) 8.197: 3 (0.000 / 0.999) 8.365: 5 (0.001 / 1.000) 8.532: 4 (0.000 / 1.000) var : p Probabilities: false: 0.9970000000000000 true: 0.0030000000000000 mean = 0.003 Percentiles: (0.001 false) (0.01 false) (0.025 false) (0.05 false) (0.1 false) (0.25 false) (0.5 false) (0.75 false) (0.84 false) (0.9 false) (0.95 false) (0.975 false) (0.99 false) (0.999 true) (0.9999 true) (0.99999 true) Histogram (total 10000) false : 9970 ############################################################ (0.997 / 0.997) true : 30 (0.003 / 1.000) */ go2 ?=> reset_store, run_model(10_000,$model2,[show_probs_trunc,mean, show_percentiles, show_histogram ]), nl, % show_store_lengths,nl, % fail, nl. go2 => true. model2() => N = 10, % Get the _coefficients* for each value 1..N TTT = [ 1/N : I in 1..N], Sample = multinomial_dist(N,TTT), % Multiply each coefficient with its value Mults = [I*Sample[I] : I in 1..N], M = Mults.mean, P = check(M >= 8), add("m",M), add("p",P). /* Get 100 values from Model 2 and try to estimate the parameters for a gaussian distribution. Num accepted samples: 1000 Total samples: 16721 (0.060%) var : mu Probabilities (truncated): 6.575317730463723: 0.0010000000000000 6.49922353052498: 0.0010000000000000 6.492756570918838: 0.0010000000000000 6.486236749443336: 0.0010000000000000 ......... 4.524865906930001: 0.0010000000000000 4.500853128033156: 0.0010000000000000 4.480986452885432: 0.0010000000000000 4.479614507630288: 0.0010000000000000 mean = 5.5065 Percentiles: (0.001 4.4809850809401768) (0.01 4.5931041249507594) (0.025 4.6849665143224248) (0.05 4.7273137558844471) (0.1 4.8128236377671474) (0.25 5.082040735325795) (0.5 5.4860002503199503) (0.75 5.964694843145411) (0.84 6.1292778848341092) (0.9 6.231447611111891) (0.95 6.3173850992775451) (0.975 6.3873009429719767) (0.99 6.444213788744162) (0.999 6.4992996247249168) (0.9999 6.5677159198898476) (0.99999 6.5745575494063431) var : sigma Probabilities (truncated): 2.092240393204726: 0.0010000000000000 1.940198550904262: 0.0010000000000000 1.898466184222356: 0.0010000000000000 1.892999451091979: 0.0010000000000000 ......... 0.005291148556998: 0.0010000000000000 0.00498036854201: 0.0010000000000000 0.004828339444859: 0.0010000000000000 0.000383153092295: 0.0010000000000000 mean = 0.886156 Percentiles: (0.001 0.0048238942585065473) (0.01 0.017433627470132721) (0.025 0.040748489108285164) (0.05 0.080567100821326995) (0.1 0.14927530621610366) (0.25 0.41724007130471996) (0.5 0.9075692707707963) (0.75 1.3475035003374813) (0.84 1.4990346165835087) (0.9 1.5932751500947751) (0.95 1.7071019756640782) (0.975 1.7712983162590759) (0.99 1.8353164857185056) (0.999 1.940350592746559) (0.9999 2.0770514131589199) (0.99999 2.0907214952001598) var : post Probabilities (truncated): 10.366263629124362: 0.0010000000000000 10.00522832651469: 0.0010000000000000 9.280754610406387: 0.0010000000000000 9.057675193274113: 0.0010000000000000 ......... 1.429896422463819: 0.0010000000000000 1.014644629949184: 0.0010000000000000 0.516968599808548: 0.0010000000000000 0.415491239497569: 0.0010000000000000 mean = 5.47227 Percentiles: (0.001 0.51686712244823707) (0.01 2.4729877670371385) (0.025 3.1968756726837246) (0.05 3.6064867615634419) (0.1 4.1285927761161973) (0.25 4.8285385519576902) (0.5 5.4902897501395493) (0.75 6.1488644490483573) (0.84 6.4276192143238067) (0.9 6.8433416291753604) (0.95 7.2686427710812671) (0.975 7.7609463990117762) (0.99 8.3378763630767256) (0.999 10.005589361817291) (0.9999 10.33019620239368) (0.99999 10.362656886451326) */ go3 ?=> run_model(100,model2,[presentation=[]]), Data = get_store().get("m"), println(data=Data), println([min=Data.min,mean=Data.mean,max=Data.max,stdev=Data.stdev]), reset_store, run_model(10_000,$model3(Data),[show_probs_trunc,mean, show_percentiles, min_accepted_samples=1000,show_accepted_samples=false ]), nl, % show_store_lengths,nl, % fail, nl. go3 => true. model3(Data) => Mu = uniform(0,10), Sigma = uniform(0,5), X = normal_dist_n(Mu,Sigma,Data.len), observe_abc(Data,X), Post = normal_dist(Mu,Sigma), if observed_ok then add("mu",Mu), add("sigma",Sigma), add("post",Post), end.