From 1e786170793775025dcbe7c7757e21475bdd89cf Mon Sep 17 00:00:00 2001 From: NunoSempere Date: Sun, 21 Jan 2024 11:44:28 +0100 Subject: [PATCH] streamline fermi estimation a little --- .../core/06_dissolving_fermi_paradox/example | Bin 21792 -> 21752 bytes .../06_dissolving_fermi_paradox/example.c | 130 +++++------------- .../core/06_dissolving_fermi_paradox/naive.c | 79 +++++++++++ 3 files changed, 111 insertions(+), 98 deletions(-) create mode 100644 examples/core/06_dissolving_fermi_paradox/naive.c diff --git a/examples/core/06_dissolving_fermi_paradox/example b/examples/core/06_dissolving_fermi_paradox/example index 45bb810266f8df999fb6c0dd08d9f906a32afc6c..21a9fdc6353b71620bd75b25080e7c6a4e3c0999 100755 GIT binary patch delta 2910 zcmZ8j3s6+o89w*I$}Y%0V3&Qdu*h8%F|wc}TBXBdg^N$57@<(3#?bLGjiwlNVhRc| zVZ;&rXBtg286a7+elBNZ2Lu@G9NPof9(JKU&d^Xh-Bx| z9>u3iUwB07(WR5z^#Ri}x>U0Za!GR^*jli&9=3sR&?Zph4GqkU7wN)18ltV!9!r{hjM;T{m(8mcEV;4cU=1r&e7AHuYvAR zw=FZmn39Bkqiv0;433oF!k7;p=YzxAU}pc|8K~GY+q4|`#+GqEpF%FjlT#dg0?CK$ z9E$p&th!1UF?OdD!%gVO2TdC|ESP9=Vc#-%|D~QLi#lM-Y~ju>7zq6ZntoxyBIrB;ts&Qv06#8L(`T_mQT`fGIPu)wzy*>(fUKg)< z=>t!;<-V?c&Pz8u*$LfVXb0|-!~2x*Ep7Mi5~Cg(@@88mbhSesP3<`k)NV~siO(lY zP=ham_PB3W;-H%!+uGzoH+|u&&)upG>psR=ER*cmja@JDn=&TbX7yl-=at{w)7>t0 zQ@vkGcntc=$$PZl|5M=}UGZnT?sB<(nU*W@5|XV~`Qk2l)3aRr(y&AM{aUi(!0VDX zfuQYaPI@|hfymP6jdaPPnb(>A>7u`;OVj@3g8TV87m++L!nH=O#7oGaVEtgxU2+6g z3gfgsknQQmn_%15fp&6Aj9@5RF$iD^rb-L z%$Sz*JUBgoGNeoJ2xv>;*X9yLeTbZVFO{Ck7@n4ZClAzDsL5sMW1*j+kJZ)TCe@jK zo=VSVdVJsT5d1ZPBTQXctr7OXA5!ReW`;PPLcN)*GD}mqGrr@9(5oy(s%n9`RUni? z`N6V;nb1+59-|k6xsEYa9cFc9X>sU7pn; z3M2G&_6cE_TBLT7k|`;)Qn*F+p$4=qQW1@Z3PU-Jl{{J5{$rJ!ryP@dH>F-+J44TCwVRKT|bI)#AzhTqXO-;MW5%wnqXNKl%Z48BS z5*h6d&oG|ma1V8br;ByP)E};CK8xFe(`Ut(Rdp0N0QBG;83l%cEdy8qu=A>_eg%2U zHC5fk<-iUgyRNDyf!ly*fnEQB9VN61u zhbmK6!Hnjfe<3iZje?|NO{izUcg9HIY}5|aFcxe^9JbhM-!GP1J*{z-)^%|wghZ>S z%xW*On#+e?z!3#kn)RbLfHPsyOt#XEf@I%KbgQsYTY$W9JqGzZ*T_0I zZAv=2J^$w2^7-g~jx8+It;6E=L@i(t1~uR&u&q0|46|z$E3Dx-2PVd6wIqqJfFDVX zLi^{=YUZcWHl(VHLFeAV~mJ9Rl1t$WrFMvnX(t6sKUSI|I@nR3e%~wna-<;IC_W3kXr~AKwZT`VQTO6J!AY0U21zVG zl%J_7YB|h%GImDukfrmpG3I5B!3QNeU!G^d!_P}dpp=T~Ik7)L?P7rDJI$ShgItNR z*lE2=>nq&eHx21JyjZ+w7XWyYJ}q|(=e71T>aB1GE*qKyc+Gi;#5smF)H|Go(SoO3O`1>$%UwErxWRGukb`CPr5* ztrnwiEWZ{fq|&}sc@|^mk4p5%Rpr6}^;WsbQ|)oYQY40FsiZnDr7NGu$4h0PB1ZRB ceP4xA5-+OCQo0;C3 zch9}Q^SkHVd(OFcVe%lm{}3CnnFI@?v^-TXHLqq2L|N6z0V|1R=SGcAQKzH95;l`f z5wg_jqq1LxnSQd5J8Ry$ZqI?vn2qTpTN6&nmKj^pa^}7k{GraQb1<6qmgVXKgSIYY z2_A|#8{JXQkNnf?KkYp(t4-(i7c%t$?{58I62o}?BOrus^QLNXVZE1cIwNIAFpTtb zkU{zDL&_gP@rW}iB_#2b^1Mp<^&=EPr3~poStH9G)J#C-C?;%*7l2&=Rjw}*BYutk zAQ0)p2a4y$N6fla>g&`jt9oAb2*s#Uwo%Gu)yjgT4&{uBXG7TYZD-iCmvV#P-sX)P zMBz4HxdCV7XP&{_BihiHiEeDo!owG42ogB`_)g;Uh(@J1pzo{aGxfbN%M&hPG12Nn$!4DK1%OK-xtnBtD?s)v?Dm+2Z7Gr zNj>nHll-oC5He+*5tTX&?DVr86`(j5s-~+@8r=%S%M5e_v1dn7@NIM|3PLarCF zMtb>c#MwV0Wm4JI@68HTJ$$jImN@I`)h5Jz&-5!sDspYG zAX|DCGt77Z)8Ap>-J?H7#iGeJ3xe82$1jE2qT~TSxs&=mrAhGp3O{G(V)98`0#~@0 z+>LY2Dg?;1&o1|N<)+f1*MFmJns;|@IwsOb_hebx0H+`5or_7HS@f2B{rA0HL&zyI zV?c5CNvA`H^K6dzYNIEam!a#A$tlbJ@rVt1D_!Wy|s8 zrqvaBEgo#=Br(OWXTKy{<9BdhSjnw;SJWMqcGjxUj^SmLhVPJkn}fYWs%;+hZ`o!= zX8=bnOwX{AOSY=nhw)@%D}TaSmx)G%b;k?Jlj!AVg@B)wem0?+-EARNb}_OMeCeGD zvdg}fogi23F3UL88<#0{s6?hux-?GgjugvDGh&i;P7NcBO+N^LY&Dxntz!|FXC`kt z#3*HXr9Qz-{_GG_@5jOXLKRb*KG7xY4_w@xZo z0HbrlNZxWLaXBW^=Ukb*6Sb~J93J|lC-C$sF{G=Ubw-kvR1mod9I1Ac{2*zzX%u(o zMTLYc8Kjp-$#0XMb8g@jE7L?No8D#rNI23)eF-02p4`o`!{pHO%aaiMAq8B?qmTB9t8ZZg* zd$c`#9{K}JV%6t;j&@UG|#clfCzUF%W@s{k4Cpe(jd@~%rQDA8zWRAayTWv8tC30A)n1FilE8;@fgXTzr;ZEc@+61V|j7Gkl1P<#)2u+ zKo4nu|FR^dES>~%Q`q$2vE0llY&yA?Ki3}qBNxnvO~;rSiDWDp$u89vhgEcM98S{t z2tHfGLl&~Mpon9vq^DpVuJ~BNc23A5TRg>Zcg&N?I>|N9_7aC`FNdBEEder?GckM9 z)p@w@rZ~Kf;*Dy)yitgmNK3dsCyiQ|V%c2`BI%a3kYjjHN0b z?%&h+HF&9MQ;Hw*_vCz0Jo_Q-b8diR;_Uykc1&A4pgbrH6Rv9p7f2wzTpa5lW5w|~ zb1gD$vs)dQAoA(dqdqK?v{M|+j*D_g+klk;UW{tYZ%LAC$@ 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)); + } + */ +}