squiggle.c/examples/more/06_nuclear_recovery/example.c

89 lines
3.9 KiB
C
Raw Normal View History

#include "../../../squiggle.h"
#include "../../../squiggle_more.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
2023-08-11 12:26:29 +00:00
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.
2023-08-11 12:26:29 +00:00
}
double yearly_probability_nuclear_collapse_2023(uint64_t* seed)
{
return yearly_probability_nuclear_collapse(2023, seed);
2023-08-11 12:26:29 +00:00
}
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 - 1);
return sample_laplace(successes, failures, seed);
2023-08-11 12:26:29 +00:00
}
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;
2023-08-11 12:26:29 +00:00
}
int main()
{
// set randomness seed
uint64_t* seed = malloc(sizeof(uint64_t));
*seed = 1000; // xorshift can't start with 0
2023-08-11 12:26:29 +00:00
int num_samples = 1000000;
// Before a first nuclear collapse
printf("## Before the first nuclear collapse\n");
ci ci_90_2023 = sampler_get_90_ci(yearly_probability_nuclear_collapse_2023, 1000000, seed);
2023-09-23 22:08:51 +00:00
printf("90%% confidence interval: [%f, %f]\n", ci_90_2023.low, ci_90_2023.high);
2023-08-11 12:26:29 +00:00
2023-12-09 19:00:43 +00:00
double* yearly_probability_nuclear_collapse_2023_samples = malloc(sizeof(double) * (size_t)num_samples);
2023-08-11 12:26:29 +00:00
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));
2023-08-11 12:26:29 +00:00
// After the first nuclear collapse
printf("\n## After the first nuclear collapse\n");
ci ci_90_2070 = sampler_get_90_ci(yearly_probability_nuclear_collapse_after_recovery_example, 1000000, seed);
2023-09-23 22:08:51 +00:00
printf("90%% confidence interval: [%f, %f]\n", ci_90_2070.low, ci_90_2070.high);
2023-08-11 12:26:29 +00:00
2023-12-09 19:00:43 +00:00
double* yearly_probability_nuclear_collapse_after_recovery_samples = malloc(sizeof(double) * (size_t)num_samples);
2023-08-11 12:26:29 +00:00
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));
2023-08-11 12:26:29 +00:00
// After the first nuclear collapse (antiinductive)
printf("\n## After the first nuclear collapse (antiinductive)\n");
ci ci_90_antiinductive = sampler_get_90_ci(yearly_probability_nuclear_collapse_after_recovery_antiinductive, 1000000, seed);
2023-09-23 22:08:51 +00:00
printf("90%% confidence interval: [%f, %f]\n", ci_90_antiinductive.low, ci_90_antiinductive.high);
2023-08-11 12:26:29 +00:00
2023-12-09 19:00:43 +00:00
double* yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples = malloc(sizeof(double) * (size_t)num_samples);
2023-08-11 12:26:29 +00:00
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(yearly_probability_nuclear_collapse_2023_samples);
free(yearly_probability_nuclear_collapse_after_recovery_samples);
free(yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples);
free(seed);
}