diff --git a/examples/08_nuclear_war/example b/examples/08_nuclear_war/example index d5a3462..4787456 100755 Binary files a/examples/08_nuclear_war/example and b/examples/08_nuclear_war/example differ diff --git a/examples/08_nuclear_war/example.c b/examples/08_nuclear_war/example.c index 785a37d..d68233f 100644 --- a/examples/08_nuclear_war/example.c +++ b/examples/08_nuclear_war/example.c @@ -1,13 +1,15 @@ #include "../../squiggle.h" #include +#include #include #include double probability_of_dying_nuno(uint64_t* seed) { double first_year_russian_nuclear_weapons = 1953; - double current_year = 2023; - double laplace_probability_nuclear_exchange_next_year = sample_beta(current_year - first_year_russian_nuclear_weapons, 0, seed); + double current_year = 2022; + double laplace_probability_nuclear_exchange_year = sample_beta(1, current_year - first_year_russian_nuclear_weapons + 1, seed); + double laplace_probability_nuclear_exchange_month = 1 - pow(1-laplace_probability_nuclear_exchange_year,(1.0/12.0)) ; double london_hit_conditional_on_russia_nuclear_weapon_usage = sample_beta(7.67, 69.65, seed); // I.e., a beta distribution with a range of 0.05 to 0.16 into: https://nunosempere.com/blog/2023/03/15/fit-beta/ @@ -17,7 +19,7 @@ double probability_of_dying_nuno(uint64_t* seed) // 0.2 to 0.8, i.e., 20% to 80%, again using the previous tool double proportion_which_die_if_bomb_drops_in_london = sample_beta(10.00, 2.45, seed); // 60% to 95% - double probability_of_dying = laplace_probability_nuclear_exchange_next_year * london_hit_conditional_on_russia_nuclear_weapon_usage * informed_actor_not_able_to_escape * proportion_which_die_if_bomb_drops_in_london; + 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; return probability_of_dying; } @@ -32,11 +34,36 @@ double probability_of_dying_eli(uint64_t* seed) return probability_of_dying; } +double mixture(uint64_t* seed){ + double (*samplers[])(uint64_t*) = {probability_of_dying_nuno, probability_of_dying_eli}; + double weights[] = {0.5, 0.5}; + return sample_mixture(samplers, weights, 2, seed); +} + int main() { // set randomness seed uint64_t* seed = malloc(sizeof(uint64_t)); *seed = 1000; // xorshift can't start with 0 + int n = 1000 * 1000; + + + double* mixture_result = malloc(sizeof(double) * n); + for(int i=0; i