diff --git a/examples/core/06_dissolving_fermi_paradox/example b/examples/core/06_dissolving_fermi_paradox/example index 45bb810..21a9fdc 100755 Binary files a/examples/core/06_dissolving_fermi_paradox/example and b/examples/core/06_dissolving_fermi_paradox/example differ diff --git a/examples/core/06_dissolving_fermi_paradox/example.c b/examples/core/06_dissolving_fermi_paradox/example.c index 8252569..ee092f8 100644 --- a/examples/core/06_dissolving_fermi_paradox/example.c +++ b/examples/core/06_dissolving_fermi_paradox/example.c @@ -4,88 +4,23 @@ #include #include -#define VERBOSE 0 - -double sample_loguniform(double a, double b, uint64_t* seed) -{ - return exp(sample_uniform(log(a), log(b), seed)); -} - int main() { // Replicate , and in particular the red line in page 11. - // Could also be interesting to just produce and save many samples. // set randomness seed uint64_t* seed = malloc(sizeof(uint64_t)); - *seed = UINT64_MAX / 64; // xorshift can't start with a seed of 0 + *seed = 1001; // xorshift can't start with a seed of 0 - // Do this naïvely, without worrying that much about numerical precision - double sample_fermi_naive(uint64_t * seed) - { - double rate_of_star_formation = sample_loguniform(1, 100, seed); - double fraction_of_stars_with_planets = sample_loguniform(0.1, 1, seed); - double number_of_habitable_planets_per_star_system = sample_loguniform(0.1, 1, seed); - double rate_of_life_formation_in_habitable_planets = sample_lognormal(1, 50, seed); - double fraction_of_habitable_planets_in_which_any_life_appears = -expm1(-rate_of_life_formation_in_habitable_planets); - // double fraction_of_habitable_planets_in_which_any_life_appears = 1-exp(-rate_of_life_formation_in_habitable_planets); - // but with more precision - double fraction_of_planets_with_life_in_which_intelligent_life_appears = sample_loguniform(0.001, 1, seed); - double fraction_of_intelligent_planets_which_are_detectable_as_such = sample_loguniform(0.01, 1, seed); - double longevity_of_detectable_civilizations = sample_loguniform(100, 10000000000, seed); - - if(VERBOSE) printf(" rate_of_star_formation = %lf\n", rate_of_star_formation); - if(VERBOSE) printf(" fraction_of_stars_with_planets = %lf\n", fraction_of_stars_with_planets); - if(VERBOSE) printf(" number_of_habitable_planets_per_star_system = %lf\n", number_of_habitable_planets_per_star_system); - if(VERBOSE) printf(" rate_of_life_formation_in_habitable_planets = %.16lf\n", rate_of_life_formation_in_habitable_planets); - if(VERBOSE) printf(" fraction_of_habitable_planets_in_which_any_life_appears = %lf\n", fraction_of_habitable_planets_in_which_any_life_appears); - if(VERBOSE) printf(" fraction_of_planets_with_life_in_which_intelligent_life_appears = %lf\n", fraction_of_planets_with_life_in_which_intelligent_life_appears); - if(VERBOSE) printf(" fraction_of_intelligent_planets_which_are_detectable_as_such = %lf\n", fraction_of_intelligent_planets_which_are_detectable_as_such); - if(VERBOSE) printf(" longevity_of_detectable_civilizations = %lf\n", longevity_of_detectable_civilizations); - - // Expected number of civilizations in the Milky way; - // see footnote 3 (p. 5) - double n = rate_of_star_formation * fraction_of_stars_with_planets * number_of_habitable_planets_per_star_system * fraction_of_habitable_planets_in_which_any_life_appears * fraction_of_planets_with_life_in_which_intelligent_life_appears * fraction_of_intelligent_planets_which_are_detectable_as_such * longevity_of_detectable_civilizations; - - return n; - } - - double sample_are_we_alone_naive(uint64_t * seed) - { - double n = sample_fermi_naive(seed); - return ((n > 1) ? 1 : 0); - } - - double n = 1000000; - double naive_fermi_proportion = 0; - for (int i = 0; i < n; i++) { - double result = sample_are_we_alone_naive(seed); - if(VERBOSE) printf("result: %lf\n", result); - naive_fermi_proportion += result; - } - printf("Naïve %% that we are not alone: %lf\n", naive_fermi_proportion / n); - - // Thinking in log space - // Taking care of numerical precision double sample_fermi_logspace(uint64_t * seed) { + // You can see a simple version of this function in naive.c in this same folder double log_rate_of_star_formation = sample_uniform(log(1), log(100), seed); double log_fraction_of_stars_with_planets = sample_uniform(log(0.1), log(1), seed); double log_number_of_habitable_planets_per_star_system = sample_uniform(log(0.1), log(1), seed); - double log_fraction_of_planets_with_life_in_which_intelligent_life_appears = sample_uniform(log(0.001), log(1), seed); - double log_fraction_of_intelligent_planets_which_are_detectable_as_such = sample_uniform(log(0.01), log(1), seed); - double log_longevity_of_detectable_civilizations = sample_uniform(log(100), log(10000000000), seed); - - if(VERBOSE) printf(" log_rate_of_star_formation = %lf\n", log_rate_of_star_formation); - if(VERBOSE) printf(" log_fraction_of_stars_with_planets = %lf\n", log_fraction_of_stars_with_planets); - if(VERBOSE) printf(" log_number_of_habitable_planets_per_star_system = %lf\n", log_number_of_habitable_planets_per_star_system); - if(VERBOSE) printf(" log_fraction_of_planets_with_life_in_which_intelligent_life_appears = %lf\n", log_fraction_of_planets_with_life_in_which_intelligent_life_appears); - if(VERBOSE) printf(" log_fraction_of_intelligent_planets_which_are_detectable_as_such = %lf\n", log_fraction_of_intelligent_planets_which_are_detectable_as_such); - if(VERBOSE) printf(" log_longevity_of_detectable_civilizations = %lf\n", log_longevity_of_detectable_civilizations); - - double log_n1 = log_rate_of_star_formation + log_fraction_of_stars_with_planets + log_number_of_habitable_planets_per_star_system + log_fraction_of_planets_with_life_in_which_intelligent_life_appears + log_fraction_of_intelligent_planets_which_are_detectable_as_such + log_longevity_of_detectable_civilizations; - if(VERBOSE) printf("first part of calculation: %lf\n", log_n1); + double log_rate_of_life_formation_in_habitable_planets = sample_normal(1, 50, seed); + double log_fraction_of_habitable_planets_in_which_any_life_appears; /* Consider: a = underlying normal @@ -93,15 +28,12 @@ int main() c = 1 - exp(-b) = fraction_of_habitable_planets_in_which_any_life_appears d = log(c) - Now, is there some way we can get d more efficiently/precisely? - - Turns out there is! - Looking at the Taylor expansion for c = 1 - exp(-b), it's b - b^2/2 + b^3/6 - x^b/24, etc. + Looking at the Taylor expansion for c = 1 - exp(-b), it's + b - b^2/2 + b^3/6 - x^b/24, etc. When b ~ 0 (as is often the case), this is close to b. - But now, if b ~ 0 - c ~ b + But now, if b ~ 0, c ~ b and d = log(c) ~ log(b) = log(exp(a)) = a Now, we could play around with estimating errors, @@ -109,15 +41,20 @@ int main() we could compute this as e.g., a < (nlog(10) + log(2))/2 so for example if we want ten digits of precision, that's a < -11 - But more empirically, the two numbers do become really close around 11 or so, and at 38 that calculation results in a -inf (so probably an overflow.) + Empirically, the two numbers as calculated in C do become really close around 11 or so, + and at 38 that calculation results in a -inf (so probably a floating point error or similar.) So we should be using that formula for somewhere between -38 << a < -11 - I chose -16 for the sake of it after playing with: - - */ - double log_rate_of_life_formation_in_habitable_planets = sample_normal(1, 50, seed); - if(VERBOSE) printf("log_rate_of_life_formation_in_habitable_planets: %lf\n", log_rate_of_life_formation_in_habitable_planets); - double log_fraction_of_habitable_planets_in_which_any_life_appears; + I chose -16 as a happy medium after playing around with + double invert(double x){ + return log(1-exp(-exp(-x))); + } + for(int i=0; i<64; i++){ + double j = i; + printf("for %lf, log(1-exp(-exp(-x))) is calculated as... %lf\n", j, invert(j)); + } + and + */ if (log_rate_of_life_formation_in_habitable_planets < -16) { log_fraction_of_habitable_planets_in_which_any_life_appears = log_rate_of_life_formation_in_habitable_planets; } else { @@ -125,36 +62,33 @@ int main() double fraction_of_habitable_planets_in_which_any_life_appears = -expm1(-rate_of_life_formation_in_habitable_planets); log_fraction_of_habitable_planets_in_which_any_life_appears = log(fraction_of_habitable_planets_in_which_any_life_appears); } - if(VERBOSE) printf(" log_fraction_of_habitable_planets_in_which_any_life_appears: %lf\n", log_fraction_of_habitable_planets_in_which_any_life_appears); - double log_n = log_n1 + log_fraction_of_habitable_planets_in_which_any_life_appears; + double log_fraction_of_planets_with_life_in_which_intelligent_life_appears = sample_uniform(log(0.001), log(1), seed); + double log_fraction_of_intelligent_planets_which_are_detectable_as_such = sample_uniform(log(0.01), log(1), seed); + double log_longevity_of_detectable_civilizations = sample_uniform(log(100), log(10000000000), seed); + + double log_n = log_rate_of_star_formation + log_fraction_of_stars_with_planets + log_number_of_habitable_planets_per_star_system + log_fraction_of_habitable_planets_in_which_any_life_appears + log_fraction_of_planets_with_life_in_which_intelligent_life_appears + log_fraction_of_intelligent_planets_which_are_detectable_as_such + log_longevity_of_detectable_civilizations; return log_n; } double sample_are_we_alone_logspace(uint64_t * seed) { - double n = sample_fermi_logspace(seed); - return ((n > 0) ? 1 : 0); + double log_n = sample_fermi_logspace(seed); + return ((log_n > 0) ? 1 : 0); + // log_n > 0 => n > 1 } double logspace_fermi_proportion = 0; - for (int i = 0; i < n; i++) { + int n_samples = 1000 * 1000; + for (int i = 0; i < n_samples; i++) { double result = sample_are_we_alone_logspace(seed); - if(VERBOSE) printf("result: %lf\n", result); logspace_fermi_proportion += result; } - printf("Using more accurate logspace computations, %% that we are not alone: %lf\n", logspace_fermi_proportion / n); + double p_not_alone = logspace_fermi_proportion / n_samples; + printf("Probability that we are not alone: %lf (%.lf%%)\n", p_not_alone, p_not_alone * 100); free(seed); - /* - double invert(double x){ - return log(1-exp(-exp(-x))); - } - for(int i=0; i<64; i++){ - double j = i; - printf("for %lf, log(1-exp(-exp(-x))) is calculated as... %lf\n", j, invert(j)); - } - */ } + diff --git a/examples/core/06_dissolving_fermi_paradox/naive.c b/examples/core/06_dissolving_fermi_paradox/naive.c new file mode 100644 index 0000000..0990a43 --- /dev/null +++ b/examples/core/06_dissolving_fermi_paradox/naive.c @@ -0,0 +1,79 @@ +#include "../../../squiggle.h" +#include +#include +#include +#include + +#define VERBOSE 0 + +double sample_loguniform(double a, double b, uint64_t* seed) +{ + return exp(sample_uniform(log(a), log(b), seed)); +} + +int main() +{ + // Replicate , and in particular the red line in page 11. + // Could also be interesting to just produce and save many samples. + + // set randomness seed + uint64_t* seed = malloc(sizeof(uint64_t)); + *seed = UINT64_MAX / 64; // xorshift can't start with a seed of 0 + + // Do this naïvely, without worrying that much about numerical precision + double sample_fermi_naive(uint64_t * seed) + { + double rate_of_star_formation = sample_loguniform(1, 100, seed); + double fraction_of_stars_with_planets = sample_loguniform(0.1, 1, seed); + double number_of_habitable_planets_per_star_system = sample_loguniform(0.1, 1, seed); + double rate_of_life_formation_in_habitable_planets = sample_lognormal(1, 50, seed); + double fraction_of_habitable_planets_in_which_any_life_appears = -expm1(-rate_of_life_formation_in_habitable_planets); + // double fraction_of_habitable_planets_in_which_any_life_appears = 1-exp(-rate_of_life_formation_in_habitable_planets); + // but with more precision + double fraction_of_planets_with_life_in_which_intelligent_life_appears = sample_loguniform(0.001, 1, seed); + double fraction_of_intelligent_planets_which_are_detectable_as_such = sample_loguniform(0.01, 1, seed); + double longevity_of_detectable_civilizations = sample_loguniform(100, 10000000000, seed); + + if(VERBOSE) printf(" rate_of_star_formation = %lf\n", rate_of_star_formation); + if(VERBOSE) printf(" fraction_of_stars_with_planets = %lf\n", fraction_of_stars_with_planets); + if(VERBOSE) printf(" number_of_habitable_planets_per_star_system = %lf\n", number_of_habitable_planets_per_star_system); + if(VERBOSE) printf(" rate_of_life_formation_in_habitable_planets = %.16lf\n", rate_of_life_formation_in_habitable_planets); + if(VERBOSE) printf(" fraction_of_habitable_planets_in_which_any_life_appears = %lf\n", fraction_of_habitable_planets_in_which_any_life_appears); + if(VERBOSE) printf(" fraction_of_planets_with_life_in_which_intelligent_life_appears = %lf\n", fraction_of_planets_with_life_in_which_intelligent_life_appears); + if(VERBOSE) printf(" fraction_of_intelligent_planets_which_are_detectable_as_such = %lf\n", fraction_of_intelligent_planets_which_are_detectable_as_such); + if(VERBOSE) printf(" longevity_of_detectable_civilizations = %lf\n", longevity_of_detectable_civilizations); + + // Expected number of civilizations in the Milky way; + // see footnote 3 (p. 5) + double n = rate_of_star_formation * fraction_of_stars_with_planets * number_of_habitable_planets_per_star_system * fraction_of_habitable_planets_in_which_any_life_appears * fraction_of_planets_with_life_in_which_intelligent_life_appears * fraction_of_intelligent_planets_which_are_detectable_as_such * longevity_of_detectable_civilizations; + + return n; + } + + double sample_are_we_alone_naive(uint64_t * seed) + { + double n = sample_fermi_naive(seed); + return ((n > 1) ? 1 : 0); + } + + double n = 1000000; + double naive_fermi_proportion = 0; + for (int i = 0; i < n; i++) { + double result = sample_are_we_alone_naive(seed); + if(VERBOSE) printf("result: %lf\n", result); + naive_fermi_proportion += result; + } + printf("Naïve %% that we are not alone: %lf\n", naive_fermi_proportion / n); + + free(seed); + + /* + double invert(double x){ + return log(1-exp(-exp(-x))); + } + for(int i=0; i<64; i++){ + double j = i; + printf("for %lf, log(1-exp(-exp(-x))) is calculated as... %lf\n", j, invert(j)); + } + */ +}