From bb91d78859dc446063287f7dee48f7771ce86499 Mon Sep 17 00:00:00 2001 From: NunoSempere Date: Sat, 20 Jan 2024 14:02:59 +0100 Subject: [PATCH] add fermi paradox example --- scratchpad/scratchpad | Bin 26912 -> 26992 bytes scratchpad/scratchpad.c | 139 +++++++++++++++++++++++++++++++++++++--- 2 files changed, 131 insertions(+), 8 deletions(-) diff --git a/scratchpad/scratchpad b/scratchpad/scratchpad index 0faebb6d265d1aa598bc8e5d14da98384e2c933d..5817d418da3e840a5704f4845ded4ab9897244ae 100755 GIT binary patch delta 6177 zcmZ`-3s@9amaeL%p+TUVM?V0OrfIMt*fa`CjDABaLK^`=K|~Z$14`lp9egBX0_vF9 zNth5%CTn))+i2ouG^-P`#BbFklSY#)$<8EbVtir5N5(Ru64a=pncDx>Bg4A8_xq~O zx&Qg=+*{|~d+x2;ev}_L%D3C1Jh7Zy@)f~#yqx}rFngF%vxmKW${n^ZO{wimqiEjk zT^FiUHOZNjIu8+y%T=P(UGx5@6a}R|KpS{#RG!MpQ(1YmoBuMo&#Q1Al*P)3*AUn3u z0o_Rc0G-oW`OEVCO^OQ2qP|5zR1lQOZ>5z%S^RNkq5L2l^#s}Yf0KFj^8MfC`)b+l z(c5Sjt~$Ms-y-AP)TLjquV2@+ESn|-f5V4RR7h6KHMGoP=YRQC#4U*Xo7@xDw;=8b z)yAM#WP32DN2eJ-wy`mSCGnr*2hV!Vg z*9<+ebTnjDcC^Z4;?8`f!UmSbO)F(89s+=;c?K)JXb6v}8-nVn4d@3zn}(pl^dZoG zpshnt16>6=JTyF_eFz#(iJ?|12D*I+8c8KU*8}Ysf|}^*P-|S2);-9;5NizW#WTHl zre_Fh?mQQIMZeYLFvL&ixIrm3p)b)6O^nF_lr86v~DNC^m&6FiMF+ zX;3J#2b0_iTvr2BsB87mj8*y-K*W}V2no!SA8_kZxO?zc;%$ZEdVpe3D1JuyTA@TJ z6to=V)EfRQQ)_rWW4;aMpwTqfX&W?F&vmxgeQmiljAC?k121rzaIbK>K4LQ+-86p= z<9=2F|1hpYxC)#G+=HCBb;*s1v8)4a6Gwq8Zuv=)HeKQ2C&+2s{Jd<*O@KxADf<#H zf`f4oO%s^1>EpcYK+U|v#w2cQjAhm=2L`Dv;%qV#6Af@6nKw@3xSQ>$?uO!P)4##Q z!*`m^2wpAIHwT(Sr9NX)E4s?TlPtl1Qs(Ebnnp`w5`^4oW_IuxRpWau}nCjul>jvreZc>lrVLPn}4L;amnUYsR&SZj}lY3i>L-bd~7V;h~xiBQ{XE+_s-JK{@)4{y;A@hP^R z6eO07t=Pf_jdLR$cotpJJDil6kS1(Rp^AiL;SVXaHerPEIZ(x|HBBYP?_HJQ2+Cx_U; z5p{z!Y(z%%+CfRW{A-zou`-B@2dQYplZj5d>?YhPj>aU!wA7+)(udkR&y09M5DKo- z35D&;YL4QanuNDR{>6>%_%<9>-~~c>_1=ZzvlD*r?|48&ZHbq zZedHO&AcO?(Puq$*p^@1?=GVMxB5BpUd>$E9%OdOp~OA{_wNjK6byk@1Ij8XK|$VJiPY;3Q1 z#=IC3#8F8$x;`>l_}W6zu1w)u94B1K!X^v&@$|Io+4ROsk~C%W#gA|Z;Jn7Ez1SYa zN?>v?r}{tf`p<5qRe!M@XRGDOTU`@`gd3z#_QM!FidvV+)vo#2ffCOKeIMF z%v76^8nggsJ+iq{X-|eju*9OtDq-J6N`L58;p#=|dZ<=tyGY5IqsN`RAdg_hKNi9v z@e>6LXW8QjrYNRUZA1WP1&leHH|9Oj^n9jc?7A2RVa{N}&t(Z0t5LkCJPZ}xsbD6S zoABi4a(?Q)K<6{(#^%D?ftyH0e_}VE|E#$FJk1zYXgQ9JT|9eU8>(B-#HRD~$5FPh z5ir50xi*@5Mr{x-m}sV(gi;fU?)6zc*vZ#K!<_pCv<73=4dZjo$k5*X6Ilao5yy>~ z3`dLomoOj3n255f1y3Jc%SsG887cdXh%=Gd^B=;JNczIFfBGMRDh~HRgTq1Nb8pCd z#$0C;lCEKiAC8Sl`P@NBNO_KcF(tJJ&nUTYGsnFt!V@ z!j3xLs^YeZSCBKkv!?pL@m5at-||*?{r4;3vKR2o;&is_7{u34cAm}I5xO-$sAX0C zFIO*Zwl}O=vCQ63ziL_I{N)Q8*VK=*_O(E@4x#*biP@ZLO|{qH#J|81HafVz4Y4iI}!qlJvvbH_8W;A{}Mb-B`7r9 zTQ!YZ^IdOte5A}R`XclQ%60%{S0$+zvi-UwnS$X6NxK!Y z4{{G=8McBGkO9d1td3nL51#`8$X4Pfk4)@Ctg=|R3qsh5nsegH!L|R-l4PeTlT%tc z{sQ03q!#uq+Kg)#!x2@j42SX}uDTd zKM)SA1H=Br5WXE)IEFq|;d_K53ZMC!0L+f*@PFB6!Dz+%S_5WE1r}z#+6DF=!+~`GW07TrVIE+;z@C&5 z<<$y+`H&kPfjdr@$Vb|%5jkME>R?-XhGtDE4Bi8^>@2+j^$ygX$c9i0FvIqrqw$4~ z;H6NvpQnXTUxB*i0zF??*wP2}18hP#vB`yv@d3h0alv{L0u8iI-x5n2E z{mNHc&~Z`RuCa8!+|Bn<)C_B!XR>M;B6OUYbCGuj4lXk+IxCk)D`!~I6aS2tM%S_7 z8d4IV3z(7Z58!SUN3CV>vLugMvrJA*%YWTBIc)9uAGm8XIDhCh$jK@ebCp~FOeU6Gx z({Ps1P_)%FPsQyla47Q1fn2ZR0enicBCj0C->G<=Rfc0YIR`=Ctr|R8bZVAc=rdDT zRiO|_rWsXkAwQiqR?VT0t72Q8!rMlni(J0!fWvZ<3;TX&9p{pjdX9>}s^KXLkcoC! z`BxP$Q;ACynA@k~b!x}87wS>)XwIW9oqF|wKI*Enwy=!Cnr~=fQ?oiaSyyl zp^E3LancfTIdIk?Q=Q^RyU$bkb?T~=H!&C&%DBd&gF;y1dX-v+uWnYJR73p&!|7ci zCh_+wV|%W`*rEt-m&V`^(vfOw+@~sUK>e|&ZEVMgpJUZu*VrXI5}zD(yW&V2lVK{} zm!#lR)OM^Y-juBX+HcvBz#VcAC#wvPU15}~PtIUCV^Cm@Eo`CY$J}{-&T3t>{QopR zzo>ZXq{;I~yGOfoX;$rvQNLWk&WrPxu2{T!d1J$hRm(ch)gBZ&w=byUc^C1EvZ4my zh!stykr${%O{|)9(}qP^p*|Hwo?SF9v|Fup)7?eme$L2RX5L2L$69zhbv>3vo_f>P zdQVgtK1W#5q=~%sS<%|CYI2H`HrBhtwEY1&_h5Y?axUMH6{Te7xw| zS)tmWT^Ai_7)Koq9gvYR;XbB delta 5655 zcmZ`-3sh9sxjyH>3^Nm$2QbV40s{;PP_E9)ucpPiMO3SE z>JGzi3hfi3>AJ3-YV@Dd*sz>75h#mS=cZQ>t59xO9-Yl6lxt9B40>6X>p4lqqMMr| zI8MizFmYzTBCZ+R5W;ExT4PyX;G)$)X~8!zXK4M(Hb-$2)p%=yZBaB7_C&5tQImgf zEmAscGg+*(Hdn>N0Pyr%gM)S(!XuhSpdoYy=xLy^Nd*3|f{Q z19TP;<9o!@2>Ou4R(Pv)(!h%2rJdK_@=vvJq7M6TkZtXLfh&M{SKefph!%|UCYuX@ zl>l;kWF_(D_mYGm#6^sDylu$ssaX)1A6RgtaU=PnV;<6{Fq4Q4$&3^@;YG6lu9V}3 zH=_GFm|&OgI3ts8wc0Viml*#-PpeVtSJ!bwTTfTPE?GZsM=8A)9WOML(t+rN%m#qm ziZPYJj)JL-4R&5^4G*p>C`LzYzAH&`gx=8=ly0`WsntGCsLY~WcDLZoqL1v42uH@y zPhyfxhwn(zkfmC;V+_^Dq=vmP2G;#|=%tu);kS3_e9QwvUnXTcUK7S;(*HS9Oy2Jl zeM=_Ybfkt2Wk7$%ca-d0C`5cmYn-cukPN11l>V+Lo}<4x++l}3&`Y>Y<*~8)v9~1& zVe&ht2xtCI`(iUqlOP?kER2}sp^LFKWBxRnb*taQHkG=q2+OKx8ZFm)hJu#N9^i+? z|GlO7eswfG8s`ps2KwU3Tl7j?wXpRTU5`tNG0Ga&Qnvw4LfwrP&mr-zHz^}N*?v47 z!OHCgUc{u;5eER!8hbQ789zxVOQ${Ysip^SDkF?Zr&ICnlnAJZO$t`zcEU=wog3=L z=DTVxKXa3kN6i!xZqnnUQk@%7&MnWiqgD=C=12Sjy&;QL_lq>zJt|$8okpLGN*2P> z=;o+I%U5nV7XGFTcEL^Vgp`bTfGH!0@rtC#ZH0$Yw-aGF4$-k8%iM@Rx~V1M+mxfJ zvi3=~5(1BIw(8zZ?q1+pCkSs0($^_(2yb1Zjj5f&Kmx_PX9{)S(xdJgel2ZC z{VD0vf>!COB$akvDDC_<*m+buf0Y)beU{A1c(BWE9L3Hf4yF%`7zIt`Y-U`g%Jej0 z_7z&4?h{;BXis{LrSOU>E`gtzE94xVAv}783P%qJzxkRxo($8eIC-D)UX7y~PrC4W zEOmH-LUk;aWo6TCPhI4H#j-l~to97{du6*{3GM|U&U zX1ps((zF-8Jcx^j3s_1{v#RxbVDcj~J{!~NptEDrLX%^VvRg9JC@U*PIAO<-RYKY&+LiT+Fyj&hvKIYSrqny6QeuACn#47`IJf}Ej3IyVm zFSV_z1Sj?nP>$DWI*WS-OSHceMP=SKLSZC*>?L83jhcL`y%Ql^83}RT1aOTO?+wd~ zt*E_${i$q#u81kPUvp8;zd}B|YNM#!YGK=XT9%t&`W0?-7)GqI(JyoVC>*oU%(44t z#Q;?tMuUUHAou(Ezwj}wvH$j!fU>&(VrCkw%h2u)QbdezP`o!jQ zR5|YVlWU-PIRT|!zVd$bRLgMPKjY6b{3-?&kDaCB@ymp{XXy_h;u%WMdpK;iC zUOsl3P3UbCnV{4(4NoW57hs;B)tu~4bdD1 z9R!UI;o7(|gnT^SNqRtopxjwWnh)9zx}4$XBxyb9;CV?p3A%AWlH8~WgGl}5pxh-% z+5&3D!S4axcv+Idb+7|86H=gl$`rnr%D?L`2 z${(Xmg^sqK4=SmFJ%wgh z!0eoAwoWx0*dTA?aba&BW_D;!JSXw22Bu?U1kEmCpUG^PexLk0KFo2mxaHWyM!z$f%622AD#){jm8S290j!ynU2MKSvMV4a`R{}d(a+rhS;pbv{O^?Sj- z@i~286sKp&^U3FwI4y4BAlT?HB&mUokH*#!idOZtaviGs)JZRvnxfx$UIIGiUK0+BY*P;ufaD z8eVo=IJ(~CAfKs-_fueY6!i$9R5q)<`{XPqk72HtIbt<+#Yp)L;f$PzlFJ?ZPW;+K z?&F%MynH(EqkZUEI6tT(&Oh=iOf?VTY@8>LE|ll+J!G5hh&2?Bbga3zb3V++f%OJv z+d~|1_&9og>YTk{tWOP?W)}}eKoZy9gXgF?K4u@7las4`<7#G4t72^*G<+=>9-p{A z+BIjous4;4=H&3>DY?Rn*_2l}VzoafG=pK8VFqg|92wf56q+;+hlj_>*Y-)nkJ6!v zB*BGHR20W?Yz=7aq_FlMFp9*I_Ku^el`94RIC`_PIL^;%1X>L(!r* zxom;kr3&_D(-#lr2(@@Mt11@!E}C7HBlL`>=c?w>a8-2MO01Smusun3AlGlkb!d39 zg0nir+GZ8+9jo-46vFLL@qTqd+Qk2?;(Ny{iVdp%2Qtpl(JEWpc~!^iP=wkVE}?7X z;y7omYNV1ORiDF$HEXCOa!b~5H3{XriAhjaA{2Zf*M85WlU0KrH7Hs^ z&H~ORWUW;kX}=RyeZP9xMd7xoYU^Ly_7r8p?go$o^e{$v#M**0)I)xYZH|Lrm5o$sraB8wUg=$ z1)cJ5SyVzII5~`}Ww~mZ&2Y9#8l2}e6{uQ+QQeU>D?;c$7q08hUvxk~S<`h7F0SQy zE0x!IZ6~3}nnqo;38+@+G@3%hJ5ljSb2^%HFzVnebOjCP>(m{n6z?_+GuJlMr|x_^xCvqqR}GS1XQb- z8ud|MqnG?keB@hVHELyUB9$!3H)^#Cm2*k{4|j!sma_gH>RTf5d#SG}-#9o=!OF?j PoNo-&D!n*FCh7kH769Ih diff --git a/scratchpad/scratchpad.c b/scratchpad/scratchpad.c index 64932b6..bd763c1 100644 --- a/scratchpad/scratchpad.c +++ b/scratchpad/scratchpad.c @@ -5,21 +5,144 @@ #include #include +double sample_loguniform(double a, double b, uint64_t* seed){ + return exp(sample_uniform(log(a), log(b), seed)); + +} + int main() { // set randomness seed uint64_t* seed = malloc(sizeof(uint64_t)); - *seed = 1000; // xorshift can't start with a seed of 0 + *seed = UINT64_MAX/64; // xorshift can't start with a seed of 0 - int n = 1000000; - double* xs = malloc(sizeof(double) * n); - for (int i = 0; i < n; i++) { - xs[i] = sample_to(10, 100, seed); + double 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); + + // printf(" rate_of_star_formation = %lf\n", rate_of_star_formation); + // printf(" fraction_of_stars_with_planets = %lf\n", fraction_of_stars_with_planets); + // printf(" number_of_habitable_planets_per_star_system = %lf\n", number_of_habitable_planets_per_star_system); + // printf(" rate_of_life_formation_in_habitable_planets = %.16lf\n", rate_of_life_formation_in_habitable_planets); + // printf(" fraction_of_habitable_planets_in_which_any_life_appears = %lf\n", fraction_of_habitable_planets_in_which_any_life_appears); + // printf(" fraction_of_planets_with_life_in_which_intelligent_life_appears = %lf\n", fraction_of_planets_with_life_in_which_intelligent_life_appears); + // printf(" fraction_of_intelligent_planets_which_are_detectable_as_such = %lf\n", fraction_of_intelligent_planets_which_are_detectable_as_such); + // 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; } - ci ci_90 = array_get_90_ci(xs, n); - printf("Recovering confidence interval of sample_to(10, 100):\n low: %f, high: %f\n", ci_90.low, ci_90.high); + double fermi_paradox_naive(uint64_t* seed){ + double n = fermi_naive(seed); + return (n > 1 ? 1 : 0); + } + + double result; + for(int i=0; i<1000; i++){ + result = fermi_naive(seed); + printf("result from fermi_naive: %lf\n", result); + printf("\n\n"); + } + printf("result from naïve implementation: %lf\n", result); + + // Thinking in log space + double fermi_logspace(uint64_t* seed){ + 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); + + // printf(" log_rate_of_star_formation = %lf\n", log_rate_of_star_formation); + // printf(" log_fraction_of_stars_with_planets = %lf\n", log_fraction_of_stars_with_planets); + // printf(" log_number_of_habitable_planets_per_star_system = %lf\n", log_number_of_habitable_planets_per_star_system); + // 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); + // printf(" log_fraction_of_intelligent_planets_which_are_detectable_as_such = %lf\n", log_fraction_of_intelligent_planets_which_are_detectable_as_such); + // 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; + printf("first part of calculation: %lf\n", log_n1); + + /* Consider fraction_of_habitable_planets_in_which_any_life_appears separately. + Imprecisely, we could do: + + double rate_of_life_formation_in_habitable_planets = sample_lognormal(1, 50, seed); + double fraction_of_habitable_planets_in_which_any_life_appears = 1- exp(-rate_of_life_formation_in_habitable_planets); + double log_fraction_of_habitable_planets_in_which_any_life_appears = log(1-fraction_of_habitable_planets_in_which_any_life_appears); + double n = exp(log_n1) * fraction_of_habitable_planets_in_which_any_life_appears; + // or: + double n2 = exp(log_n1 + log(fraction_of_habitable_planets_in_which_any_life_appears)) + + However, we lose all precision here. + + Now, say + a = underlying normal + b = rate_of_life_formation_in_habitable_planets = exp(underlying normal) + c = 1 - exp(-b) = fraction_of_habitable_planets_in_which_any_life_appears + d = log(c) + + Now, is there some way we can 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. + When b ~ 0 (as it is), this is close to b. + + But now, if b ~ 0 + c ~ b + and d = log(c) ~ log(b) = log(exp(a)) = a + */ + double log_rate_of_life_formation_in_habitable_planets = sample_normal(1, 50, seed); + 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; + if(log_rate_of_life_formation_in_habitable_planets < -32){ + log_fraction_of_habitable_planets_in_which_any_life_appears = log_rate_of_life_formation_in_habitable_planets; + } else{ + double rate_of_life_formation_in_habitable_planets = exp(log_rate_of_life_formation_in_habitable_planets); + 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); + } + 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; + + return log_n; + } + + double result2; + + /* + for(int i=0; i<1000; i++){ + result2 = fermi_logspace(seed); + printf("result from logspace implementation: %lf.2\n", result2); + printf("\n\n"); + } + */ - printf("Size of uint64_t: %ld", sizeof(uint64_t*)); free(seed); }