/* Hypothesis testing in Picat. From the AgenaRisk model Tutorial/Hypothesis Testing Comparison of two materials A and B which has different number of tests of faultiness: - A was tested in 200 cases where 10 was faulty - B was tested in 100 cases where 9 was fault. Is A better then B? (This is - yet another variant of an - A/B test.) Cf my Gamble model gamble_hypothesis_testing.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. main => go. /* Note: This does not give many samples (see the lengths below). var : a is faulty Probabilities (truncated): 0.094123993912927: 0.0086206896551724 0.093227403638998: 0.0086206896551724 0.086424585557063: 0.0086206896551724 0.08445250967704: 0.0086206896551724 ......... 0.024231908889275: 0.0086206896551724 0.023184212734994: 0.0086206896551724 0.020687234745118: 0.0086206896551724 0.01794079800619: 0.0086206896551724 mean = 0.0510554 HPD intervals: HPD interval (0.5): 0.03088996329867..0.05032713245489 HPD interval (0.84): 0.03246421309565..0.07359039481503 HPD interval (0.9): 0.02465432896548..0.07359039481503 HPD interval (0.95): 0.02068723474512..0.08163495577496 HPD interval (0.99): 0.02068723474512..0.09412399391293 HPD interval (0.999): 0.01794079800619..0.09412399391293 HPD interval (0.9999): 0.01794079800619..0.09412399391293 HPD interval (0.99999): 0.01794079800619..0.09412399391293 Histogram: 0.018: 1 ##### (0.009 / 0.009) 0.020: 1 ##### (0.009 / 0.017) 0.022: 0 (0.000 / 0.017) 0.024: 2 ########### (0.017 / 0.034) 0.026: 2 ########### (0.017 / 0.052) 0.027: 2 ########### (0.017 / 0.069) 0.029: 1 ##### (0.009 / 0.078) 0.031: 2 ########### (0.017 / 0.095) 0.033: 8 ############################################ (0.069 / 0.164) 0.035: 3 ################ (0.026 / 0.190) 0.037: 6 ################################# (0.052 / 0.241) 0.039: 3 ################ (0.026 / 0.267) 0.041: 3 ################ (0.026 / 0.293) 0.043: 5 ########################### (0.043 / 0.336) 0.045: 2 ########### (0.017 / 0.353) 0.047: 11 ############################################################ (0.095 / 0.448) 0.048: 7 ###################################### (0.060 / 0.509) 0.050: 8 ############################################ (0.069 / 0.578) 0.052: 0 (0.000 / 0.578) 0.054: 5 ########################### (0.043 / 0.621) 0.056: 3 ################ (0.026 / 0.647) 0.058: 3 ################ (0.026 / 0.672) 0.060: 5 ########################### (0.043 / 0.716) 0.062: 5 ########################### (0.043 / 0.759) 0.064: 2 ########### (0.017 / 0.776) 0.066: 8 ############################################ (0.069 / 0.845) 0.067: 1 ##### (0.009 / 0.853) 0.069: 1 ##### (0.009 / 0.862) 0.071: 5 ########################### (0.043 / 0.905) 0.073: 4 ###################### (0.034 / 0.940) 0.075: 0 (0.000 / 0.940) 0.077: 0 (0.000 / 0.940) 0.079: 2 ########### (0.017 / 0.957) 0.081: 1 ##### (0.009 / 0.966) 0.083: 0 (0.000 / 0.966) 0.085: 1 ##### (0.009 / 0.974) 0.087: 1 ##### (0.009 / 0.983) 0.088: 0 (0.000 / 0.983) 0.090: 0 (0.000 / 0.983) 0.092: 2 ########### (0.017 / 1.000) var : b is faulty Probabilities (truncated): 0.15257823296988: 0.0086206896551724 0.146706141768399: 0.0086206896551724 0.146681110784707: 0.0086206896551724 0.146370362818936: 0.0086206896551724 ......... 0.046607086152435: 0.0086206896551724 0.046596186414887: 0.0086206896551724 0.045142276506392: 0.0086206896551724 0.041092760715708: 0.0086206896551724 mean = 0.0876653 HPD intervals: HPD interval (0.5): 0.07252497708688..0.10084014933302 HPD interval (0.84): 0.04514227650639..0.10938748710694 HPD interval (0.9): 0.04978100036937..0.12615460113893 HPD interval (0.95): 0.04109276071571..0.13162616335042 HPD interval (0.99): 0.04109276071571..0.14670614176840 HPD interval (0.999): 0.04109276071571..0.15257823296988 HPD interval (0.9999): 0.04109276071571..0.15257823296988 HPD interval (0.99999): 0.04109276071571..0.15257823296988 Histogram: 0.041: 1 ##### (0.009 / 0.009) 0.044: 1 ##### (0.009 / 0.017) 0.047: 2 ########## (0.017 / 0.034) 0.049: 2 ########## (0.017 / 0.052) 0.052: 3 ############### (0.026 / 0.078) 0.055: 1 ##### (0.009 / 0.086) 0.058: 4 #################### (0.034 / 0.121) 0.061: 3 ############### (0.026 / 0.147) 0.063: 3 ############### (0.026 / 0.172) 0.066: 3 ############### (0.026 / 0.198) 0.069: 3 ############### (0.026 / 0.224) 0.072: 7 ################################### (0.060 / 0.284) 0.075: 7 ################################### (0.060 / 0.345) 0.077: 4 #################### (0.034 / 0.379) 0.080: 5 ######################### (0.043 / 0.422) 0.083: 5 ######################### (0.043 / 0.466) 0.086: 6 ############################## (0.052 / 0.517) 0.088: 6 ############################## (0.052 / 0.569) 0.091: 3 ############### (0.026 / 0.595) 0.094: 5 ######################### (0.043 / 0.638) 0.097: 3 ############### (0.026 / 0.664) 0.100: 12 ############################################################ (0.103 / 0.767) 0.102: 2 ########## (0.017 / 0.784) 0.105: 2 ########## (0.017 / 0.802) 0.108: 5 ######################### (0.043 / 0.845) 0.111: 1 ##### (0.009 / 0.853) 0.114: 1 ##### (0.009 / 0.862) 0.116: 3 ############### (0.026 / 0.888) 0.119: 2 ########## (0.017 / 0.905) 0.122: 1 ##### (0.009 / 0.914) 0.125: 2 ########## (0.017 / 0.931) 0.127: 1 ##### (0.009 / 0.940) 0.130: 2 ########## (0.017 / 0.957) 0.133: 0 (0.000 / 0.957) 0.136: 1 ##### (0.009 / 0.966) 0.139: 0 (0.000 / 0.966) 0.141: 0 (0.000 / 0.966) 0.144: 0 (0.000 / 0.966) 0.147: 3 ############### (0.026 / 0.991) 0.150: 1 ##### (0.009 / 1.000) var : hypothesis Probabilities: A better than B: 0.8965517241379310 A not better than B: 0.1034482758620690 mean = [A better than B = 0.896552,A not better than B = 0.103448] HPD intervals: show_hpd_intervals: data is not numeric Histogram: A better than B: 104 ############################################################ (0.897 / 0.897) A not better than B: 12 ####### (0.103 / 1.000) var : p(a is faulty) < p(b is faulty) Probabilities: true: 0.8965517241379310 false: 0.1034482758620690 mean = [true = 0.896552,false = 0.103448] HPD intervals: show_hpd_intervals: data is not numeric Histogram: false : 12 ####### (0.103 / 0.103) true : 104 ############################################################ (0.897 / 1.000) Variable lengths: a is faulty = 116 b is faulty = 116 hypothesis = 116 p(a is faulty) < p(b is faulty) = 116 */ go ?=> reset_store, run_model(100_000,$model,[show_probs_trunc,mean,show_hpd_intervals,show_histogram]), nl, show_store_lengths(), % fail, nl. go => true. model() => ATests = 200, BTests = 100, ACount = 10, BCount = 9, ProbAIsFaulty = beta_dist(1,ACount), ProbBIsFaulty = beta_dist(1,BCount), AFaults = binomial_dist(ATests,ProbAIsFaulty), % AFaults = bern_sum(ProbAIsFaulty,ATests), BFaults = binomial_dist(BTests,ProbBIsFaulty), % BFaults = bern_sum(ProbBIsFaulty,BTests), Hypothesis = cond(ProbAIsFaulty < ProbBIsFaulty, "A better than B", "A not better than B"), observe(AFaults==ACount,BFaults=BCount), if observed_ok then add("a is faulty",ProbAIsFaulty), add("b is faulty",ProbBIsFaulty), add("p(a is faulty) < p(b is faulty)",cond(ProbAIsFaulty