example.c (2833B)
1 #include "../../../squiggle.h" 2 #include "../../../squiggle_more.h" 3 #include <math.h> 4 #include <stdio.h> 5 #include <stdlib.h> 6 7 double probability_of_dying_nuno(uint64_t* seed) 8 { 9 double first_year_russian_nuclear_weapons = 1953; 10 double current_year = 2022; 11 double laplace_probability_nuclear_exchange_year = sample_beta(1, current_year - first_year_russian_nuclear_weapons + 1, seed); 12 double laplace_probability_nuclear_exchange_month = 1 - pow(1 - laplace_probability_nuclear_exchange_year, (1.0 / 12.0)); 13 14 double london_hit_conditional_on_russia_nuclear_weapon_usage = sample_beta(7.67, 69.65, seed); 15 // I.e., a beta distribution with a range of 0.05 to 0.16 into: https://nunosempere.com/blog/2023/03/15/fit-beta/ 16 // 0.05 were my estimate and Samotsvety's estimate in March 2022, respectively: 17 // https://forum.effectivealtruism.org/posts/KRFXjCqqfGQAYirm5/samotsvety-nuclear-risk-forecasts-march-2022#Nu_o_Sempere 18 double informed_actor_not_able_to_escape = sample_beta(3.26212166586967, 3.26228162008564, seed); 19 // 0.2 to 0.8, i.e., 20% to 80%, again using the previous tool 20 double proportion_which_die_if_bomb_drops_in_london = sample_beta(10.00, 2.45, seed); // 60% to 95% 21 22 double probability_of_dying = laplace_probability_nuclear_exchange_month * london_hit_conditional_on_russia_nuclear_weapon_usage * informed_actor_not_able_to_escape * proportion_which_die_if_bomb_drops_in_london; 23 return probability_of_dying; 24 } 25 26 double probability_of_dying_eli(uint64_t* seed) 27 { 28 double russia_nato_nuclear_exchange_in_next_month = sample_beta(1.30, 1182.99, seed); // .0001 to .003 29 double london_hit_conditional = sample_beta(3.47, 8.97, seed); // 0.1 to 0.5 30 double informed_actors_not_able_to_escape = sample_beta(2.73, 5.67, seed); // .1 to .6 31 double proportion_which_die_if_bomb_drops_in_london = sample_beta(3.00, 1.46, seed); // 0.3 to 0.95; 32 33 double probability_of_dying = russia_nato_nuclear_exchange_in_next_month * london_hit_conditional * informed_actors_not_able_to_escape * proportion_which_die_if_bomb_drops_in_london; 34 return probability_of_dying; 35 } 36 37 double sample_nuclear_model(uint64_t* seed) 38 { 39 double (*samplers[])(uint64_t*) = { probability_of_dying_nuno, probability_of_dying_eli }; 40 double weights[] = { 0.5, 0.5 }; 41 return sample_mixture(samplers, weights, 2, seed); 42 } 43 44 int main() 45 { 46 // set randomness seed 47 uint64_t* seed = malloc(sizeof(uint64_t)); 48 *seed = 1000; // xorshift can't start with 0 49 50 int n = 1 * MILLION; 51 double* xs = malloc(sizeof(double) * (size_t)n); 52 for (int i = 0; i < n; i++) { 53 xs[i] = sample_nuclear_model(seed); 54 } 55 56 printf("\n# Stats\n"); 57 array_print_stats(xs, n); 58 printf("\n# Histogram\n"); 59 array_print_90_ci_histogram(xs, n, 20); 60 61 free(xs); 62 free(seed); 63 }