diff --git a/examples/08_nuclear_war/example b/examples/08_nuclear_war/example index 4787456..8df8418 100755 Binary files a/examples/08_nuclear_war/example and b/examples/08_nuclear_war/example differ diff --git a/examples/10_nuclear_recovery/example b/examples/10_nuclear_recovery/example new file mode 100755 index 0000000..d33589a Binary files /dev/null and b/examples/10_nuclear_recovery/example differ diff --git a/examples/10_nuclear_recovery/example.c b/examples/10_nuclear_recovery/example.c index 88a4746..670355d 100644 --- a/examples/10_nuclear_recovery/example.c +++ b/examples/10_nuclear_recovery/example.c @@ -8,13 +8,37 @@ double laplace(double successes, double trials){ return (successes + 1)/(trials + 2); } -double yearly_probability_nuclear_apocalypse(double year, uint64_t* seed) +double yearly_probability_nuclear_collapse(double year, uint64_t* seed) { double successes = 0; double failures = (year - 1960); - + return sample_laplace(successes, failures, seed); + // ^ can change to (successes + 1)/(trials + 2) + // to get a probability, + // rather than sampling from a distribution over probabilities. +} +double yearly_probability_nuclear_collapse_2023(uint64_t* seed){ + return yearly_probability_nuclear_collapse(2023, seed); } +double yearly_probability_nuclear_collapse_after_recovery(double year, double rebuilding_period_length_years, uint64_t* seed){ + // assumption: nuclear + double successes = 1.0; + double failures = (year - rebuilding_period_length_years - 1960); + return sample_laplace(successes, failures, seed); +} +double yearly_probability_nuclear_collapse_after_recovery_example(uint64_t* seed){ + double year = 2070; + double rebuilding_period_length_years =30; + // So, there was a nuclear collapse in 2040, + // then a recovery period of 30 years + // and it's now 2070 + return yearly_probability_nuclear_collapse_after_recovery(year, rebuilding_period_length_years, seed); +} + +double yearly_probability_nuclear_collapse_after_recovery_antiinductive(uint64_t* seed){ + return yearly_probability_nuclear_collapse(2023, seed)/2; +} int main() { @@ -22,7 +46,40 @@ int main() uint64_t* seed = malloc(sizeof(uint64_t)); *seed = 1000; // xorshift can't start with 0 - int n = 1000 * 1000; + int num_samples = 1000000; + + // Before a first nuclear collapse + printf("## Before the first nuclear collapse\n"); + struct c_i c_i_90_2023 = get_90_confidence_interval(yearly_probability_nuclear_collapse_2023, seed); + printf("90%% confidence interval: [%f, %f]\n", c_i_90_2023.low, c_i_90_2023.high); + + double* yearly_probability_nuclear_collapse_2023_samples = malloc(sizeof(double) * num_samples); + for (int i = 0; i < num_samples; i++) { + yearly_probability_nuclear_collapse_2023_samples[i] = yearly_probability_nuclear_collapse_2023(seed); + } + printf("mean: %f\n", array_mean(yearly_probability_nuclear_collapse_2023_samples, num_samples)); + + // After the first nuclear collapse + printf("\n## After the first nuclear collapse\n"); + struct c_i c_i_90_2070 = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_example, seed); + printf("90%% confidence interval: [%f, %f]\n", c_i_90_2070.low, c_i_90_2070.high); + + double* yearly_probability_nuclear_collapse_after_recovery_samples = malloc(sizeof(double) * num_samples); + for (int i = 0; i < num_samples; i++) { + yearly_probability_nuclear_collapse_after_recovery_samples[i] = yearly_probability_nuclear_collapse_after_recovery_example(seed); + } + printf("mean: %f\n", array_mean(yearly_probability_nuclear_collapse_after_recovery_samples, num_samples)); + + // After the first nuclear collapse (antiinductive) + printf("\n## After the first nuclear collapse (antiinductive)\n"); + struct c_i c_i_90_antiinductive = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_antiinductive, seed); + printf("90%% confidence interval: [%f, %f]\n", c_i_90_antiinductive.low, c_i_90_antiinductive.high); + + double* yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples = malloc(sizeof(double) * num_samples); + for (int i = 0; i < num_samples; i++) { + yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples[i] = yearly_probability_nuclear_collapse_after_recovery_antiinductive(seed); + } + printf("mean: %f\n", array_mean(yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples, num_samples)); free(seed); } diff --git a/squiggle.c b/squiggle.c index 0150b6a..adb50e8 100644 --- a/squiggle.c +++ b/squiggle.c @@ -135,6 +135,7 @@ double sample_beta(double a, double b, uint64_t* seed) } double sample_laplace(double successes, double failures, uint64_t* seed){ + // see return sample_beta(successes + 1, failures + 1, seed); }