/* Wiener process in Picat. From Mathematica WienerProcess """ The state at time t follows NormalDistribution[Mu t,Sigma Sqrt[t]] """ For more, see ppl_distributions.pi and ppl_distributions_test.pi Cf my Gamble model gamble_wiener_process.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. /* Two Wiener processes, check the difference of their last (10th) values. var : v1 Probabilities (truncated): 14.705761865605325: 0.0010000000000000 13.854485128927735: 0.0010000000000000 12.678854496538444: 0.0010000000000000 12.395752064086247: 0.0010000000000000 ......... -5.944407893831976: 0.0010000000000000 -6.366442855405332: 0.0010000000000000 -7.114835067966704: 0.0010000000000000 -7.143483832676957: 0.0010000000000000 mean = 2.73941 [len = 1000,min = -7.14348,mean = 2.73941,median = 2.70832,max = 14.7058,variance = 11.7002,stdev = 3.42055] HPD intervals: HPD interval (0.84): -1.92677217350563..7.49212093026340 HPD interval (0.9): -2.55686761104214..8.41933412771617 HPD interval (0.94): -3.86084338885962..8.55286975829761 HPD interval (0.99): -5.38440486606200..12.29211582952607 HPD interval (0.99999): -7.14348383267696..14.70576186560533 var : v2 Probabilities (truncated): 14.728662002171928: 0.0010000000000000 13.229872389427939: 0.0010000000000000 12.971451633472123: 0.0010000000000000 12.897808998498629: 0.0010000000000000 ......... -6.815275028974423: 0.0010000000000000 -7.104104376904557: 0.0010000000000000 -8.052909950813191: 0.0010000000000000 -11.027484420381716: 0.0010000000000000 mean = 3.06025 [len = 1000,min = -11.0275,mean = 3.06025,median = 2.97594,max = 14.7287,variance = 12.8052,stdev = 3.57844] HPD intervals: HPD interval (0.84): -1.88433492177190..7.85804977413445 HPD interval (0.9): -2.54414521112013..8.76138902187483 HPD interval (0.94): -3.25272144875929..9.77129496704160 HPD interval (0.99): -5.53108263075066..12.97145163347212 HPD interval (0.99999): -11.02748442038172..14.72866200217193 var : diff Probabilities (truncated): 15.370884009095652: 0.0010000000000000 15.360292497726068: 0.0010000000000000 15.093684245732824: 0.0010000000000000 14.397148532426488: 0.0010000000000000 ......... -12.590202607305006: 0.0010000000000000 -13.416134297154521: 0.0010000000000000 -14.755546972109219: 0.0010000000000000 -15.666144714466: 0.0010000000000000 mean = -0.320834 [len = 1000,min = -15.6661,mean = -0.320834,median = -0.253559,max = 15.3709,variance = 25.0425,stdev = 5.00425] HPD intervals: HPD interval (0.84): -7.19382984369016..6.50797935311670 HPD interval (0.9): -7.72374498402124..8.51022423138059 HPD interval (0.94): -9.67129432299803..9.02213962079223 HPD interval (0.99): -12.59020260730501..13.45679019461571 HPD interval (0.99999): -15.66614471446600..15.37088400909565 Diff: Histogram (total 1000) -15.666: 1 # (0.001 / 0.001) -14.890: 1 # (0.001 / 0.002) -14.114: 0 (0.000 / 0.002) -13.338: 1 # (0.001 / 0.003) -12.562: 7 ###### (0.007 / 0.010) -11.787: 4 ### (0.004 / 0.014) -11.011: 9 ######## (0.009 / 0.023) -10.235: 7 ###### (0.007 / 0.030) -9.459: 12 ########## (0.012 / 0.042) -8.683: 17 ############## (0.017 / 0.059) -7.907: 10 ######## (0.010 / 0.069) -7.131: 26 ###################### (0.026 / 0.095) -6.355: 30 ######################### (0.030 / 0.125) -5.579: 38 ################################ (0.038 / 0.163) -4.803: 43 #################################### (0.043 / 0.206) -4.027: 52 ########################################### (0.052 / 0.258) -3.251: 59 ################################################# (0.059 / 0.317) -2.475: 57 ################################################ (0.057 / 0.374) -1.699: 45 ###################################### (0.045 / 0.419) -0.924: 50 ########################################## (0.050 / 0.469) -0.148: 72 ############################################################ (0.072 / 0.541) 0.628: 72 ############################################################ (0.072 / 0.613) 1.404: 46 ###################################### (0.046 / 0.659) 2.180: 58 ################################################ (0.058 / 0.717) 2.956: 60 ################################################## (0.060 / 0.777) 3.732: 51 ########################################### (0.051 / 0.828) 4.508: 26 ###################### (0.026 / 0.854) 5.284: 27 ####################### (0.027 / 0.881) 6.060: 31 ########################## (0.031 / 0.912) 6.836: 26 ###################### (0.026 / 0.938) 7.612: 18 ############### (0.018 / 0.956) 8.388: 11 ######### (0.011 / 0.967) 9.163: 7 ###### (0.007 / 0.974) 9.939: 9 ######## (0.009 / 0.983) 10.715: 3 ### (0.003 / 0.986) 11.491: 1 # (0.001 / 0.987) 12.267: 5 #### (0.005 / 0.992) 13.043: 0 (0.000 / 0.992) 13.819: 4 ### (0.004 / 0.996) 14.595: 4 ### (0.004 / 1.000) */ go ?=> reset_store, run_model(1000,$model,[show_probs_trunc,mean, show_simple_stats, % show_percentiles, show_hpd_intervals,hpd_intervals=[0.84,0.9,0.94,0.99,0.99999] % show_histogram % min_accepted_samples=1000,show_accepted_samples=true ]), println("Diff:"), Diff = get_store().get("diff"), Diff.show_histogram, nl, % show_store_lengths,nl, % fail, nl. go => true. model() => N = 10, Mu = 1, Sigma = 2, T = 3, X1 = wiener_process_dist_n(Mu,Sigma,T, N), X2 = wiener_process_dist_n(Mu,Sigma,T, N), V1 = X1.last, V2 = X2.last, Diff = V1-V2, add("v1",V1), add("v2",V2), add("diff",Diff). /* Parameter recovery [12.4476,8.07192,11.7114,11.8962,13.6799,6.05413,7.5259,11.5393,6.95803,9.47007,8.80363,10.1392,7.80246,9.11638,14.7742,13.7802,10.2021,10.9913,5.28251,9.94025,1.59735,6.58536,5.21428,11.9675,5.50749,8.13085,5.09463,4.72915,19.6298,6.3046,9.11231,12.4323,9.99145,7.14377,10.2191,9.21102,8.94428,15.6623,7.10844,11.575,9.23241,11.5585,8.12443,13.2771,7.01499,7.68817,12.0122,13.8447,3.88139,14.6755,9.22588,4.26994,5.98704,10.4594,13.2531,7.65485,3.73224,1.88725,8.63202,3.70237,6.53409,15.5478,13.4834,12.604,6.84211,9.34582,7.89809,11.4495,9.10887,10.5022,9.27473,2.9063,8.45627,14.1601,14.1333,11.0288,7.025,11.213,10.2415,7.703,10.6967,5.63914,8.47795,13.0716,12.2944,7.47932,6.42877,13.3116,7.66858,10.6851,14.9274,9.77613,10.4034,12.7218,8.47316,13.9616,8.10836,12.1025,10.6611,13.6653] [len = 100,min = 1.59735,mean = 9.56473,median = 9.31028,max = 19.6298,variance = 11.1277,stdev = 3.33582] Histogram (total 100) 1.597: 1 ####### (0.010 / 0.010) 2.048: 1 ####### (0.010 / 0.020) 2.499: 0 (0.000 / 0.020) 2.950: 1 ####### (0.010 / 0.030) 3.401: 0 (0.000 / 0.030) 3.851: 3 #################### (0.030 / 0.060) 4.302: 1 ####### (0.010 / 0.070) 4.753: 1 ####### (0.010 / 0.080) 5.204: 3 #################### (0.030 / 0.110) 5.655: 2 ############# (0.020 / 0.130) 6.105: 3 #################### (0.030 / 0.160) 6.556: 3 #################### (0.030 / 0.190) 7.007: 6 ######################################## (0.060 / 0.250) 7.458: 4 ########################### (0.040 / 0.290) 7.909: 8 ##################################################### (0.080 / 0.370) 8.360: 3 #################### (0.030 / 0.400) 8.810: 3 #################### (0.030 / 0.430) 9.261: 9 ############################################################ (0.090 / 0.520) 9.712: 1 ####### (0.010 / 0.530) 10.163: 6 ######################################## (0.060 / 0.590) 10.614: 6 ######################################## (0.060 / 0.650) 11.064: 3 #################### (0.030 / 0.680) 11.515: 5 ################################# (0.050 / 0.730) 11.966: 4 ########################### (0.040 / 0.770) 12.417: 4 ########################### (0.040 / 0.810) 12.868: 2 ############# (0.020 / 0.830) 13.318: 4 ########################### (0.040 / 0.870) 13.769: 5 ################################# (0.050 / 0.920) 14.220: 2 ############# (0.020 / 0.940) 14.671: 2 ############# (0.020 / 0.960) 15.122: 1 ####### (0.010 / 0.970) 15.572: 2 ############# (0.020 / 0.990) 16.023: 0 (0.000 / 0.990) 16.474: 0 (0.000 / 0.990) 16.925: 0 (0.000 / 0.990) 17.376: 0 (0.000 / 0.990) 17.827: 0 (0.000 / 0.990) 18.277: 0 (0.000 / 0.990) 18.728: 0 (0.000 / 0.990) 19.179: 1 ####### (0.010 / 1.000) [[mean,stdev],Data appears roughly symmetric with light tails. Using mean and stdev.] min_accepted_samples = 100 var : mu Probabilities (truncated): 12.21467913697226: 0.0100000000000000 11.784991906855717: 0.0100000000000000 11.615974461480963: 0.0100000000000000 11.535221483341987: 0.0100000000000000 ......... 1.020461228220007: 0.0100000000000000 0.975067262060506: 0.0100000000000000 0.940009579500188: 0.0100000000000000 0.800261227786663: 0.0100000000000000 mean = 5.9039 [len = 100,min = 0.800261,mean = 5.9039,median = 5.24276,max = 12.2147,variance = 10.7666,stdev = 3.28126] HPD intervals: HPD interval (0.84): 0.80026122778666..9.58030540942229 HPD interval (0.9): 0.97506726206051..10.91673849658889 HPD interval (0.94): 0.94000957950019..11.18984367288176 HPD interval (0.99): 0.80026122778666..11.78499190685572 HPD interval (0.99999): 0.80026122778666..12.21467913697226 Histogram (total 100) 0.800: 2 ################# (0.020 / 0.020) 1.086: 5 ########################################### (0.050 / 0.070) 1.371: 3 ########################## (0.030 / 0.100) 1.656: 3 ########################## (0.030 / 0.130) 1.942: 1 ######### (0.010 / 0.140) 2.227: 2 ################# (0.020 / 0.160) 2.512: 5 ########################################### (0.050 / 0.210) 2.798: 1 ######### (0.010 / 0.220) 3.083: 5 ########################################### (0.050 / 0.270) 3.369: 2 ################# (0.020 / 0.290) 3.654: 1 ######### (0.010 / 0.300) 3.939: 2 ################# (0.020 / 0.320) 4.225: 6 ################################################### (0.060 / 0.380) 4.510: 3 ########################## (0.030 / 0.410) 4.795: 5 ########################################### (0.050 / 0.460) 5.081: 3 ########################## (0.030 / 0.490) 5.366: 7 ############################################################ (0.070 / 0.560) 5.651: 1 ######### (0.010 / 0.570) 5.937: 0 (0.000 / 0.570) 6.222: 1 ######### (0.010 / 0.580) 6.507: 0 (0.000 / 0.580) 6.793: 0 (0.000 / 0.580) 7.078: 5 ########################################### (0.050 / 0.630) 7.364: 2 ################# (0.020 / 0.650) 7.649: 0 (0.000 / 0.650) 7.934: 2 ################# (0.020 / 0.670) 8.220: 3 ########################## (0.030 / 0.700) 8.505: 2 ################# (0.020 / 0.720) 8.790: 6 ################################################### (0.060 / 0.780) 9.076: 3 ########################## (0.030 / 0.810) 9.361: 2 ################# (0.020 / 0.830) 9.646: 2 ################# (0.020 / 0.850) 9.932: 2 ################# (0.020 / 0.870) 10.217: 0 (0.000 / 0.870) 10.503: 2 ################# (0.020 / 0.890) 10.788: 3 ########################## (0.030 / 0.920) 11.073: 3 ########################## (0.030 / 0.950) 11.359: 1 ######### (0.010 / 0.960) 11.644: 3 ########################## (0.030 / 0.990) 11.929: 1 ######### (0.010 / 1.000) var : sigma Probabilities (truncated): 6.867301173912036: 0.0100000000000000 5.971397797563764: 0.0100000000000000 5.5756588073334: 0.0100000000000000 5.494312795574922: 0.0100000000000000 ......... 0.115033942328316: 0.0100000000000000 0.109933326072029: 0.0100000000000000 0.011733244178692: 0.0100000000000000 0.00561529770755: 0.0100000000000000 mean = 2.34532 [len = 100,min = 0.0056153,mean = 2.34532,median = 2.13049,max = 6.8673,variance = 2.17861,stdev = 1.47601] HPD intervals: HPD interval (0.84): 0.00561529770755..3.67742139551669 HPD interval (0.9): 0.10993332607203..4.66360069562849 HPD interval (0.94): 0.00561529770755..4.87621679197821 HPD interval (0.99): 0.00561529770755..5.97139779756376 HPD interval (0.99999): 0.00561529770755..6.86730117391204 Histogram (total 100) 0.006: 2 ############ (0.020 / 0.020) 0.177: 3 ################## (0.030 / 0.050) 0.349: 4 ######################## (0.040 / 0.090) 0.520: 1 ###### (0.010 / 0.100) 0.692: 4 ######################## (0.040 / 0.140) 0.863: 2 ############ (0.020 / 0.160) 1.035: 10 ############################################################ (0.100 / 0.260) 1.206: 1 ###### (0.010 / 0.270) 1.378: 3 ################## (0.030 / 0.300) 1.549: 4 ######################## (0.040 / 0.340) 1.721: 5 ############################## (0.050 / 0.390) 1.893: 5 ############################## (0.050 / 0.440) 2.064: 7 ########################################## (0.070 / 0.510) 2.236: 4 ######################## (0.040 / 0.550) 2.407: 4 ######################## (0.040 / 0.590) 2.579: 3 ################## (0.030 / 0.620) 2.750: 3 ################## (0.030 / 0.650) 2.922: 6 #################################### (0.060 / 0.710) 3.093: 7 ########################################## (0.070 / 0.780) 3.265: 3 ################## (0.030 / 0.810) 3.436: 0 (0.000 / 0.810) 3.608: 3 ################## (0.030 / 0.840) 3.780: 0 (0.000 / 0.840) 3.951: 1 ###### (0.010 / 0.850) 4.123: 1 ###### (0.010 / 0.860) 4.294: 3 ################## (0.030 / 0.890) 4.466: 0 (0.000 / 0.890) 4.637: 3 ################## (0.030 / 0.920) 4.809: 2 ############ (0.020 / 0.940) 4.980: 1 ###### (0.010 / 0.950) 5.152: 0 (0.000 / 0.950) 5.323: 1 ###### (0.010 / 0.960) 5.495: 2 ############ (0.020 / 0.980) 5.667: 0 (0.000 / 0.980) 5.838: 0 (0.000 / 0.980) 6.010: 1 ###### (0.010 / 0.990) 6.181: 0 (0.000 / 0.990) 6.353: 0 (0.000 / 0.990) 6.524: 0 (0.000 / 0.990) 6.696: 1 ###### (0.010 / 1.000) var : t Probabilities: 1: 0.4200000000000000 2: 0.2500000000000000 3: 0.1300000000000000 4: 0.0700000000000000 9: 0.0500000000000000 7: 0.0300000000000000 10: 0.0200000000000000 6: 0.0200000000000000 5: 0.0100000000000000 mean = 2.62 [len = 100,min = 1,mean = 2.62,median = 2.0,max = 10,variance = 5.3356,stdev = 2.30989] HPD intervals: HPD interval (0.84): 1.00000000000000..4.00000000000000 HPD interval (0.9): 1.00000000000000..6.00000000000000 HPD interval (0.94): 1.00000000000000..9.00000000000000 HPD interval (0.99): 1.00000000000000..10.00000000000000 HPD interval (0.99999): 1.00000000000000..10.00000000000000 Histogram (total 100) 1: 42 ############################################################ (0.420 / 0.420) 2: 25 #################################### (0.250 / 0.670) 3: 13 ################### (0.130 / 0.800) 4: 7 ########## (0.070 / 0.870) 5: 1 # (0.010 / 0.880) 6: 2 ### (0.020 / 0.900) 7: 3 #### (0.030 / 0.930) 9: 5 ####### (0.050 / 0.980) 10: 2 ### (0.020 / 1.000) var : post Probabilities (truncated): 19.220575742013452: 0.0100000000000000 17.605816948052929: 0.0100000000000000 16.527122543306525: 0.0100000000000000 16.212570476808736: 0.0100000000000000 ......... 0.630324691215438: 0.0100000000000000 0.629974718716985: 0.0100000000000000 0.226452904715279: 0.0100000000000000 -4.252617640696839: 0.0100000000000000 mean = 9.50623 [len = 100,min = -4.25262,mean = 9.50623,median = 9.4885,max = 19.2206,variance = 13.7356,stdev = 3.70616] HPD intervals: HPD interval (0.84): 4.76334943393613..14.05012359285788 HPD interval (0.9): 4.76334943393613..15.60110331857075 HPD interval (0.94): 2.10156242943216..16.52712254330653 HPD interval (0.99): 0.22645290471528..19.22057574201345 HPD interval (0.99999): -4.25261764069684..19.22057574201345 Histogram (total 100) -4.253: 1 ##### (0.010 / 0.010) -3.666: 0 (0.000 / 0.010) -3.079: 0 (0.000 / 0.010) -2.492: 0 (0.000 / 0.010) -1.905: 0 (0.000 / 0.010) -1.318: 0 (0.000 / 0.010) -0.732: 0 (0.000 / 0.010) -0.145: 0 (0.000 / 0.010) 0.442: 3 ############## (0.030 / 0.040) 1.029: 0 (0.000 / 0.040) 1.616: 0 (0.000 / 0.040) 2.203: 2 ######### (0.020 / 0.060) 2.789: 0 (0.000 / 0.060) 3.376: 0 (0.000 / 0.060) 3.963: 0 (0.000 / 0.060) 4.550: 1 ##### (0.010 / 0.070) 5.137: 4 ################## (0.040 / 0.110) 5.723: 1 ##### (0.010 / 0.120) 6.310: 5 ####################### (0.050 / 0.170) 6.897: 5 ####################### (0.050 / 0.220) 7.484: 6 ############################ (0.060 / 0.280) 8.071: 2 ######### (0.020 / 0.300) 8.658: 13 ############################################################ (0.130 / 0.430) 9.244: 7 ################################ (0.070 / 0.500) 9.831: 5 ####################### (0.050 / 0.550) 10.418: 11 ################################################### (0.110 / 0.660) 11.005: 6 ############################ (0.060 / 0.720) 11.592: 11 ################################################### (0.110 / 0.830) 12.179: 2 ######### (0.020 / 0.850) 12.765: 1 ##### (0.010 / 0.860) 13.352: 3 ############## (0.030 / 0.890) 13.939: 2 ######### (0.020 / 0.910) 14.526: 2 ######### (0.020 / 0.930) 15.113: 1 ##### (0.010 / 0.940) 15.700: 2 ######### (0.020 / 0.960) 16.286: 2 ######### (0.020 / 0.980) 16.873: 0 (0.000 / 0.980) 17.460: 1 ##### (0.010 / 0.990) 18.047: 0 (0.000 / 0.990) 18.634: 1 ##### (0.010 / 1.000) Post: normalized_rmse = 0.195682 Excellent match (very small difference between distributions). abc_fit_score_explained = [0.237346,Good fit; only small distributional differences.] */ go2 ?=> Mu = 3, Sigma = 2, T = 3, Data = wiener_process_dist_n(Mu,Sigma,T,100), println(Data), show_simple_stats(Data), show_histogram(Data), println(recommend_abc_stats_explained(Data)), reset_store, run_model(1000,$model2(Data),[show_probs_trunc,mean, show_simple_stats, % show_percentiles, show_hpd_intervals,hpd_intervals=[0.84,0.9,0.94,0.99,0.99999], show_histogram, min_accepted_samples=100 % ,show_accepted_samples=true ]), analyse_post_simple(Data), nl, % show_store_lengths,nl, % fail, nl. model2(Data) => Mu = uniform(0,100), Sigma = uniform(0,10), T = random_integer1(10), X = wiener_process_dist_n(Mu,Sigma,T, Data.len), % observe_abc(Data,X,1/2), observe_abc(Data,X,1,[mean,stdev,q1,q25,q75,q99]), Post = wiener_process_dist(Mu,Sigma,T), if observed_ok then add("mu",Mu), add("sigma",Sigma), add("t",T), add("post",Post) end.