/* Thrice exceptional in Picat. From @cremieuxrecueil https://x.com/cremieuxrecueil/status/1858655075192713685 """ I simulated 100,000 people to show how often people are "thrice-exceptional": Smart, stable, and exceptionally hard-working. ... In short, there are not many people who are thrice-exceptional, in the sense of being at least +2 standard deviations in conscientiousness, emotional stability (i.e., inverse neuroticism), and intelligence. """ Also see Gilles E. Gignac "The number of exceptional people: Fewer than 85 per 1 million across key traits" https://www.sciencedirect.com/science/article/pii/S019188692400415X The criteria for exceptional is 2SD for each dimension. Assuming that all traits are normal distributed (with mean 0 and stddev 1), i.e. a 2SD of 2: Picat> X = 1-normal_dist_cdf(0,1,2), Inv = 1/X X = 0.022750131948179 Inv = 43.955789015985438 I.e. about one in 44. For three exceptional (2sd) traits: Picat> X = (1-normal_dist_cdf(0,1,2))**3, Inv = 1/X X = 0.00001177475175 Inv = 84927.480527094347053 About one in 850 000 people. It's very rate with three 3sd traits: Picat> X = (1-normal_dist_cdf(0,1,3))**3, Inv = 1/X X = 0.000000002459818 Inv = 406534219.626161277294159 Being thrice 1SD exceptional is quite rare (about 1 in 250) Picat> X = (1-normal_dist_cdf(0,1,1))**3, Inv = 1/X X = 0.00399358907433 Inv = 250.401326072294211 Cf my Gamble model gamble_thrice_expectional.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. /* The simulation below assumes that the traits are are all normally distributed (here with mean 0 and standard deviation 1), and does not correlate the traits (contrary to the research cited above, though the correlation seems to be small, see the discussion in the thread), It shows a little less probability of exceptional (2sd) people: about 10 in 1_000_000. Here are results from some runs: [x = [2.04341,2.14531,2.4399],p2sd = true,p3sd = false] [x = [2.19173,2.42895,2.30804],p2sd = true,p3sd = false] [x = [2.1929,2.84213,2.45515],p2sd = true,p3sd = false] [x = [2.31177,2.75989,2.15428],p2sd = true,p3sd = false] [x = [2.03581,2.39419,2.66046],p2sd = true,p3sd = false] [x = [2.61802,2.30708,2.14415],p2sd = true,p3sd = false] [x = [2.12672,2.39587,3.15556],p2sd = true,p3sd = false] [x = [2.48715,2.06973,2.16322],p2sd = true,p3sd = false] [x = [2.2924,2.45396,2.71669],p2sd = true,p3sd = false] [x = [2.28954,2.22335,2.67683],p2sd = true,p3sd = false] [x = [2.67614,2.04826,2.0434],p2sd = true,p3sd = false] Result: num_exceptional1d = 4034 num_exceptional2d = 11 num_exceptional3d = 0 Theoretical: sd = 1 = 0.00399359 = 3993.59 sd = 2 = 1.17748e-05 = 11.7748 sd = 3 = 2.45982e-09 = 0.00245982 [x = [2.44379,2.11111,2.68422],p2sd = true,p3sd = false] [x = [2.48289,2.20294,2.36165],p2sd = true,p3sd = false] [x = [2.54439,2.45125,2.05763],p2sd = true,p3sd = false] [x = [2.10683,2.6178,2.33183],p2sd = true,p3sd = false] [x = [2.15815,2.14032,2.93784],p2sd = true,p3sd = false] [x = [2.25337,2.40121,2.0901],p2sd = true,p3sd = false] [x = [2.32734,2.22259,3.14394],p2sd = true,p3sd = false] [x = [2.29623,2.82951,2.19886],p2sd = true,p3sd = false] [x = [2.01524,2.48247,2.30403],p2sd = true,p3sd = false] Result: num_exceptional1d = 4010 num_exceptional2d = 9 num_exceptional3d = 0 Theoretical: sd = 1 = 0.00399359 = 3993.59 sd = 2 = 1.17748e-05 = 11.7748 sd = 3 = 2.45982e-09 = 0.00245982 [x = [2.35901,2.37273,2.80424],p2sd = true,p3sd = false] [x = [2.17523,2.36461,2.32782],p2sd = true,p3sd = false] [x = [2.41716,2.06242,2.49106],p2sd = true,p3sd = false] [x = [2.16583,2.88526,2.09386],p2sd = true,p3sd = false] [x = [2.92118,2.48949,2.26872],p2sd = true,p3sd = false] [x = [2.66141,2.93732,2.01172],p2sd = true,p3sd = false] [x = [2.22702,2.6041,2.01904],p2sd = true,p3sd = false] [x = [2.07422,2.28967,2.52213],p2sd = true,p3sd = false] [x = [2.4409,2.06669,2.69939],p2sd = true,p3sd = false] [x = [2.14628,2.53527,2.67887],p2sd = true,p3sd = false] [x = [2.3471,2.52968,2.47936],p2sd = true,p3sd = false] [x = [3.11189,2.5253,2.49101],p2sd = true,p3sd = false] [x = [2.20969,2.20925,2.31778],p2sd = true,p3sd = false] Result: num_exceptional1d = 4094 num_exceptional2d = 13 num_exceptional3d = 0 Theoretical: sd = 1 = 0.00399359 = 3993.59 sd = 2 = 1.17748e-05 = 11.7748 sd = 3 = 2.45982e-09 = 0.00245982 And here's a run for normal_dist(100,15): [x = [133.872,134.466,132.767],p2sd = true,p3sd = false] [x = [134.764,133.518,133.957],p2sd = true,p3sd = false] [x = [134.525,134.664,132.878],p2sd = true,p3sd = false] [x = [131.713,138.045,139.922],p2sd = true,p3sd = false] [x = [132.324,131.178,131.96],p2sd = true,p3sd = false] [x = [139.801,134.329,136.74],p2sd = true,p3sd = false] [x = [135.075,134.248,136.249],p2sd = true,p3sd = false] [x = [144.149,137.094,134.419],p2sd = true,p3sd = false] [x = [137.555,130.479,136.545],p2sd = true,p3sd = false] [x = [144.635,131.883,130.995],p2sd = true,p3sd = false] [x = [132.098,152.696,136.585],p2sd = true,p3sd = false] [x = [132.645,131.82,138.768],p2sd = true,p3sd = false] [x = [130.902,133.286,138.197],p2sd = true,p3sd = false] [x = [136.814,136.279,130.651],p2sd = true,p3sd = false] [x = [130.901,141.008,136.103],p2sd = true,p3sd = false] [x = [136.578,132.581,139.105],p2sd = true,p3sd = false] Result: num_exceptional1d = 4051 num_exceptional2d = 16 num_exceptional3d = 0 Theoretical: sd = 1 = 0.00399359 = 3993.59 sd = 2 = 1.17748e-05 = 11.7748 sd = 3 = 2.45982e-09 = 0.00245982 */ go ?=> NumTraits = 3, Mu = 0, Std = 1, reset_store, Map = get_global_map(), NumRuns = 10_000, run_model(NumRuns,$model(NumTraits,Mu,Std),[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=10,show_accepted_samples=true ]), nl, println("Result:"), foreach(E in [num_exceptional1d,num_exceptional2d,num_exceptional3d]) println(E=Map.get(E,0)) end, nl, println("Theoretical:"), foreach(SD in 1..3) LimitSD = Mu + Std * SD, T = (1-normal_dist_cdf(Mu,Std,LimitSD))**NumTraits, println(sd=SD=T=T*NumRuns) end, % show_store_lengths,nl, % fail, nl. go => true. model(NumTraits,Mu,Std) => Limit1SD = Mu + Std * 1, Limit2SD = Mu + Std * 2, Limit3SD = Mu + Std * 3, X = normal_dist_n(Mu,Std,NumTraits), % Probability of n-exceptional P1SD = check( [cond(X[Trait] >= Limit1SD,1,0) : Trait in 1..NumTraits].sum == 3 ), P2SD = check( [cond(X[Trait] >= Limit2SD,1,0) : Trait in 1..NumTraits].sum == 3 ), P3SD = check( [cond(X[Trait] >= Limit3SD,1,0) : Trait in 1..NumTraits].sum == 3 ), Map = get_global_map(), if P1SD then Map.put(num_exceptional1d,Map.get(num_exceptional1d,0)+1), end, if P2SD then println([x=X,p2sd=P2SD,p3sd=P3SD]), Map.put(num_exceptional2d,Map.get(num_exceptional2d,0)+1), end, if P3SD then Map.put(num_exceptional3d,Map.get(num_exceptional3d,0)+1), end. % add("p1 sd",P1SD), % add("p2 sd",P2SD), % add("p3 sd",P3SD).