/* Student T dist in Picat. 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): 283.63676585996518: 0.0001000000000000 126.644302449597532: 0.0001000000000000 49.255110258277696: 0.0001000000000000 41.782838489617973: 0.0001000000000000 ......... -33.476611600154456: 0.0001000000000000 -38.065473958639203: 0.0001000000000000 -48.937221331682125: 0.0001000000000000 -82.178898513235737: 0.0001000000000000 [len = 10000,min = -82.1789,mean = 0.0314998,median = -0.00400102,max = 283.637,variance = 16.3302,stdev = 4.04107] Scatter plot: ymax = 283.637 | * | | | | | | | | * | | | | * | * * * |************************************************************ |************************************************************ | * * * ** * *** * * * * | ** | | * ymin = -82.1789 x = [1,10000] Histogram (total 10000) -82.179: 1 (0.000 / 0.000) -73.034: 0 (0.000 / 0.000) -63.888: 0 (0.000 / 0.000) -54.743: 0 (0.000 / 0.000) -45.597: 1 (0.000 / 0.000) -36.452: 2 (0.000 / 0.000) -27.307: 4 (0.000 / 0.001) -18.161: 10 (0.001 / 0.002) -9.016: 197 # (0.020 / 0.021) 0.130: 9585 ############################################################ (0.959 / 0.980) 9.275: 180 # (0.018 / 0.998) 18.420: 13 (0.001 / 0.999) 27.566: 3 (0.000 / 1.000) 36.711: 0 (0.000 / 1.000) 45.857: 2 (0.000 / 1.000) 55.002: 0 (0.000 / 1.000) 64.147: 0 (0.000 / 1.000) 73.293: 0 (0.000 / 1.000) 82.438: 0 (0.000 / 1.000) 91.584: 0 (0.000 / 1.000) 100.729: 0 (0.000 / 1.000) 109.874: 0 (0.000 / 1.000) 119.020: 0 (0.000 / 1.000) 128.165: 1 (0.000 / 1.000) 137.311: 0 (0.000 / 1.000) 146.456: 0 (0.000 / 1.000) 155.601: 0 (0.000 / 1.000) 164.747: 0 (0.000 / 1.000) 173.892: 0 (0.000 / 1.000) 183.037: 0 (0.000 / 1.000) 192.183: 0 (0.000 / 1.000) 201.328: 0 (0.000 / 1.000) 210.474: 0 (0.000 / 1.000) 219.619: 0 (0.000 / 1.000) 228.764: 0 (0.000 / 1.000) 237.910: 0 (0.000 / 1.000) 247.055: 0 (0.000 / 1.000) 256.201: 0 (0.000 / 1.000) 265.346: 0 (0.000 / 1.000) 274.491: 1 (0.000 / 1.000) var : d2 Probabilities (truncated): 61.387765441219877: 0.0001000000000000 37.26438696209722: 0.0001000000000000 28.628837715588993: 0.0001000000000000 27.087663136514646: 0.0001000000000000 ......... -23.347167957706294: 0.0001000000000000 -24.220566099914649: 0.0001000000000000 -25.608167193189612: 0.0001000000000000 -26.411834524255099: 0.0001000000000000 [len = 10000,min = -26.4118,mean = 0.979438,median = 0.982162,max = 61.3878,variance = 10.5566,stdev = 3.24909] Scatter plot: ymax = 61.3878 | * | | | | | * | | * * | * * * * | *** ** * * * | ** * * **** * * * * ***** ** ** | * **** ***** ******** ****** ******* ******** ***** ***** |************************************************************ |************************************************************ |************************************************************ |********* *************************** ***************** *** | * * * * ** * *** ** * ** * * ******* | * * * * * * | * * * * * ** | * * * ymin = -26.4118 x = [1,10000] Histogram (total 10000) -26.412: 2 (0.000 / 0.000) -24.217: 2 (0.000 / 0.000) -22.022: 4 (0.000 / 0.001) -19.827: 3 (0.000 / 0.001) -17.632: 1 (0.000 / 0.001) -15.437: 4 (0.000 / 0.002) -13.242: 16 (0.002 / 0.003) -11.047: 18 (0.002 / 0.005) -8.852: 52 # (0.005 / 0.010) -6.657: 139 ### (0.014 / 0.024) -4.462: 368 ####### (0.037 / 0.061) -2.267: 1184 ###################### (0.118 / 0.179) -0.072: 3295 ############################################################ (0.330 / 0.509) 2.123: 3214 ########################################################### (0.321 / 0.830) 4.318: 1118 #################### (0.112 / 0.942) 6.513: 379 ####### (0.038 / 0.980) 8.708: 106 ## (0.011 / 0.991) 10.903: 40 # (0.004 / 0.995) 13.098: 26 (0.003 / 0.997) 15.293: 11 (0.001 / 0.998) 17.488: 4 (0.000 / 0.999) 19.683: 6 (0.001 / 0.999) 21.878: 1 (0.000 / 0.999) 24.073: 1 (0.000 / 0.999) 26.268: 3 (0.000 / 1.000) 28.463: 1 (0.000 / 1.000) 30.658: 0 (0.000 / 1.000) 32.853: 0 (0.000 / 1.000) 35.048: 0 (0.000 / 1.000) 37.243: 1 (0.000 / 1.000) 39.438: 0 (0.000 / 1.000) 41.633: 0 (0.000 / 1.000) 43.828: 0 (0.000 / 1.000) 46.023: 0 (0.000 / 1.000) 48.218: 0 (0.000 / 1.000) 50.413: 0 (0.000 / 1.000) 52.608: 0 (0.000 / 1.000) 54.803: 0 (0.000 / 1.000) 56.998: 0 (0.000 / 1.000) 59.193: 1 (0.000 / 1.000) */ go ?=> reset_store, run_model(10_000,$model,[show_probs_trunc,show_simple_stats, % show_percentiles, % show_hpd_intervals,hpd_intervals=[0.94], show_scatter_plot, show_histogram % min_accepted_samples=1000,show_accepted_samples=true ]), nl, % show_store_lengths,nl, % fail, nl. go => true. model() => D1 = student_t_dist(2), D2 = student_t_dist(1,2,3), add("d1",D1), add("d2",D2). /* Parameter recovery Not very good: data = [0.375195,1.86401,-1.14067,-1.31701,0.427572,0.0397485,-2.57041,-1.65263,0.426453,-0.55268,-0.389721,0.293898,0.972209,-0.737694,-0.19484,0.70376,-0.461139,-1.71048,0.271112,-0.467443,0.43234,1.22471,4.60656,-0.667831,-0.0860493,0.584529,1.22707,-0.999334,0.933181,-0.205997,-0.138049,-1.25571,-5.18145e-05,0.0204412,-0.170974,-0.986325,-0.735796,-0.395463,-0.276406,-0.519396,-1.56763,-1.79287,4.37678,-0.0719691,1.06342,1.59937,-0.933117,0.185233,0.349776,2.09866,0.714297,-2.2762,-0.957988,-0.969876,0.620063,-2.64633,-0.179508,0.495368,0.827851,-0.3541,2.80637,0.0123938,1.2098,2.58872,1.02361,-1.30876,-1.33423,2.22744,0.570258,-2.28252,-0.351574,-1.26932,-0.719661,0.907963,1.2259,-0.111269,0.936201,-0.483348,-2.13188,-0.919186,-0.87669,-0.833298,1.39494,-2.97643,-1.05655,-0.632746,-2.24249,1.15234,0.977266,2.8561,-1.02037,-1.15742,-1.44692,1.04041,0.745445,-0.6904,1.36359,-1.66113,0.380557,-1.92032] methods = [median,iqr,mad,skewness] explain = Data is clearly skewed. Using robust statistics (median, IQR, MAD, skewness). [len = 100,min = -2.97643,mean = -0.0665531,median = -0.175241,max = 4.60656,variance = 1.89115,stdev = 1.37519] Histogram (total 100) -2.976: 1 ######## (0.010 / 0.010) -2.787: 0 (0.000 / 0.010) -2.597: 2 ############### (0.020 / 0.030) -2.408: 0 (0.000 / 0.030) -2.218: 4 ############################## (0.040 / 0.070) -2.029: 0 (0.000 / 0.070) -1.839: 2 ############### (0.020 / 0.090) -1.649: 4 ############################## (0.040 / 0.130) -1.460: 1 ######## (0.010 / 0.140) -1.270: 5 ###################################### (0.050 / 0.190) -1.081: 6 ############################################# (0.060 / 0.250) -0.891: 6 ############################################# (0.060 / 0.310) -0.702: 6 ############################################# (0.060 / 0.370) -0.512: 5 ###################################### (0.050 / 0.420) -0.322: 5 ###################################### (0.050 / 0.470) -0.133: 8 ############################################################ (0.080 / 0.550) 0.057: 4 ############################## (0.040 / 0.590) 0.246: 3 ####################### (0.030 / 0.620) 0.436: 7 ##################################################### (0.070 / 0.690) 0.625: 5 ###################################### (0.050 / 0.740) 0.815: 3 ####################### (0.030 / 0.770) 1.005: 7 ##################################################### (0.070 / 0.840) 1.194: 5 ###################################### (0.050 / 0.890) 1.384: 2 ############### (0.020 / 0.910) 1.573: 1 ######## (0.010 / 0.920) 1.763: 0 (0.000 / 0.920) 1.953: 1 ######## (0.010 / 0.930) 2.142: 2 ############### (0.020 / 0.950) 2.332: 0 (0.000 / 0.950) 2.521: 1 ######## (0.010 / 0.960) 2.711: 0 (0.000 / 0.960) 2.900: 2 ############### (0.020 / 0.980) 3.090: 0 (0.000 / 0.980) 3.280: 0 (0.000 / 0.980) 3.469: 0 (0.000 / 0.980) 3.659: 0 (0.000 / 0.980) 3.848: 0 (0.000 / 0.980) 4.038: 0 (0.000 / 0.980) 4.227: 0 (0.000 / 0.980) 4.417: 2 ############### (0.020 / 1.000) Num accepted samples: 100 Total samples: 236 (0.424%) var : eta Probabilities (truncated): 9.875923169253358: 0.0100000000000000 9.780132596278625: 0.0100000000000000 9.727412323806162: 0.0100000000000000 9.695444692250083: 0.0100000000000000 ......... 2.287591290794122: 0.0100000000000000 1.817371283572806: 0.0100000000000000 1.673036394488549: 0.0100000000000000 1.655647612482145: 0.0100000000000000 mean = 6.13521 [len = 100,min = 1.65565,mean = 6.13521,median = 6.27861,max = 9.87592,variance = 5.24464,stdev = 2.29012] var : post Probabilities (truncated): 8.171849178690854: 0.0100000000000000 3.524852989144254: 0.0100000000000000 2.793798842066005: 0.0100000000000000 2.578518557286572: 0.0100000000000000 ......... -2.540409523598509: 0.0100000000000000 -3.959811145869748: 0.0100000000000000 -4.115833123370679: 0.0100000000000000 -12.22846483589173: 0.0100000000000000 mean = 0.17987 [len = 100,min = -12.2285,mean = 0.17987,median = 0.230623,max = 8.17185,variance = 3.71449,stdev = 1.9273] Post: normalized_rmse = 0.340719 Good match; only small differences across the distribution. abc_fit_score_explained = [0.490693,Moderate mismatch; distributions differ in shape.] Mathematica: FindDistributionParameters[data, StudentTDistribution[eta]] -> {eta -> 3.68299} */ go2 ?=> Data = student_t_dist_n(3,100), println(data=Data), [Methods,Explain] = recommend_abc_stats_explained(Data), println(methods=Methods), println(explain=Explain), show_simple_stats(Data), show_histogram(Data), % show_scatter_plot(Data), reset_store, run_model(10_000,$model2(Data),[show_probs_trunc,mean,show_simple_stats, % show_percentiles, % show_hpd_intervals,hpd_intervals=[0.94], % show_scatter_plot, % show_histogram, min_accepted_samples=100,show_accepted_samples=true ]), analyse_post_simple(Data), nl, % show_store_lengths,nl, % fail, nl. go2 => true. model2(Data) => Eta = uniform(1,10), Y = student_t_dist_n(Eta,Data.len), % observe_abc(Data,Y,1/30), observe_abc(Data,Y,1/2,[median,iqr,mad,skewness]), Post = student_t_dist(Eta), if observed_ok then add("eta",Eta), add("post",Post) end.