squiggle.c

commit ef40ef5ae78a807b1c8a6e9be3a9fd18e26e4dca
parent b76a630d5893d6d64cc7fb929050c03e6cc20acc
Author: NunoSempere <nuno.sempere@protonmail.com>
Date:   Mon, 24 Jul 2023 00:14:44 +0200

change the nuclear probability to monthly

Diffstat:
M
examples/08_nuclear_war/example
 | 
0

M
examples/08_nuclear_war/example.c
 | 
33
++++++++++++++++++++++++++++++---

2 files changed, 30 insertions(+), 3 deletions(-)
diff --git a/examples/08_nuclear_war/example b/examples/08_nuclear_war/example
Binary files differ.
diff --git a/examples/08_nuclear_war/example.c b/examples/08_nuclear_war/example.c
@@ -1,13 +1,15 @@
 #include "../../squiggle.h"
 #include <stdint.h>
+#include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
 
 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<n; i++){
+      mixture_result[i] = mixture(seed);
+    }
+    
+    printf("mixture_result: [ ");
+    for(int i=0; i<9; i++){
+      printf("%.6f, ", mixture_result[i]);
+    }
+    printf("... ]\n");
+    
+		struct c_i c_i_90 = get_90_confidence_interval(mixture, seed);
+    printf("mean: %f\n", array_mean(mixture_result, n));
+    printf("90%% confidence interval: [%f, %f]\n", c_i_90.low, c_i_90.high);
+
+
     free(seed);
 }