example.c (3695B)
1 #include "../../../squiggle.h" 2 #include "../../../squiggle_more.h" 3 #include <math.h> 4 #include <stdio.h> 5 #include <stdlib.h> 6 7 double yearly_probability_nuclear_collapse(double year, uint64_t* seed) 8 { 9 double successes = 0; 10 double failures = (year - 1960); 11 return sample_laplace(successes, failures, seed); 12 // ^ can change to (successes + 1)/(trials + 2) 13 // to get a probability, 14 // rather than sampling from a distribution over probabilities. 15 } 16 double yearly_probability_nuclear_collapse_2023(uint64_t* seed) 17 { 18 return yearly_probability_nuclear_collapse(2023, seed); 19 } 20 21 double yearly_probability_nuclear_collapse_after_recovery(double year, double rebuilding_period_length_years, uint64_t* seed) 22 { 23 // assumption: nuclear 24 double successes = 1.0; 25 double failures = (year - rebuilding_period_length_years - 1960 - 1); 26 return sample_laplace(successes, failures, seed); 27 } 28 double yearly_probability_nuclear_collapse_after_recovery_example(uint64_t* seed) 29 { 30 double year = 2070; 31 double rebuilding_period_length_years = 30; 32 // So, there was a nuclear collapse in 2040, 33 // then a recovery period of 30 years 34 // and it's now 2070 35 return yearly_probability_nuclear_collapse_after_recovery(year, rebuilding_period_length_years, seed); 36 } 37 38 double yearly_probability_nuclear_collapse_after_recovery_antiinductive(uint64_t* seed) 39 { 40 return yearly_probability_nuclear_collapse(2023, seed) / 2; 41 } 42 43 int main() 44 { 45 // set randomness seed 46 uint64_t* seed = malloc(sizeof(uint64_t)); 47 *seed = 1000; // xorshift can't start with 0 48 49 int n_samples = 1000000; 50 51 // Before a first nuclear collapse 52 printf("## Before the first nuclear collapse\n"); 53 double* yearly_probability_nuclear_collapse_2023_samples = malloc(sizeof(double) * (size_t)n_samples); 54 for (int i = 0; i < n_samples; i++) { 55 yearly_probability_nuclear_collapse_2023_samples[i] = yearly_probability_nuclear_collapse_2023(seed); 56 } 57 ci ci_90_2023 = array_get_90_ci(yearly_probability_nuclear_collapse_2023_samples, n_samples); 58 printf("90%% confidence interval: [%f, %f]\n", ci_90_2023.low, ci_90_2023.high); 59 60 // After the first nuclear collapse 61 printf("\n## After the first nuclear collapse\n"); 62 63 double* yearly_probability_nuclear_collapse_after_recovery_samples = malloc(sizeof(double) * (size_t)n_samples); 64 for (int i = 0; i < n_samples; i++) { 65 yearly_probability_nuclear_collapse_after_recovery_samples[i] = yearly_probability_nuclear_collapse_after_recovery_example(seed); 66 } 67 ci ci_90_2070 = array_get_90_ci(yearly_probability_nuclear_collapse_after_recovery_samples, 1000000); 68 printf("90%% confidence interval: [%f, %f]\n", ci_90_2070.low, ci_90_2070.high); 69 70 // After the first nuclear collapse (antiinductive) 71 printf("\n## After the first nuclear collapse (antiinductive)\n"); 72 double* yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples = malloc(sizeof(double) * (size_t)n_samples); 73 for (int i = 0; i < n_samples; i++) { 74 yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples[i] = yearly_probability_nuclear_collapse_after_recovery_antiinductive(seed); 75 } 76 ci ci_90_antiinductive = array_get_90_ci(yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples, 1000000); 77 printf("90%% confidence interval: [%f, %f]\n", ci_90_antiinductive.low, ci_90_antiinductive.high); 78 79 // free seeds 80 free(yearly_probability_nuclear_collapse_2023_samples); 81 free(yearly_probability_nuclear_collapse_after_recovery_samples); 82 free(yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples); 83 free(seed); 84 }