From ef40ef5ae78a807b1c8a6e9be3a9fd18e26e4dca Mon Sep 17 00:00:00 2001 From: NunoSempere Date: Mon, 24 Jul 2023 00:14:44 +0200 Subject: [PATCH] change the nuclear probability to monthly --- examples/08_nuclear_war/example | Bin 22320 -> 22448 bytes examples/08_nuclear_war/example.c | 33 +++++++++++++++++++++++++++--- 2 files changed, 30 insertions(+), 3 deletions(-) diff --git a/examples/08_nuclear_war/example b/examples/08_nuclear_war/example index d5a34625d8c245d9bc6e48486d449e4b7d43aa40..47874565c39358521080f3b9ddceabe77de8f5fb 100755 GIT binary patch delta 5767 zcmZ`-3wRV&maeLFceU^W$%(%fUGhx?-04i}-$H4Ba4C5mpXrYnB@brPT|E&i`KlgrL z)w%aSk9+UA=bn43x?kbD|CR6XxCY`mwFEW7jdW)J*4Ay)THT$0n`$pk(;CHT6vyX= zx}vqFK{-)t^B@*H`2|`#f8h`5NkMB5(sthMl5{Camy+mZ{@n?^d78|RX$&)t62sHe z!lI1-uZw=Yi+-|2tx0&g_|$O0-Yey&#}!QeV%;-tanl78ZHvlEOn@b0v>dlxwVXf) zq9*e%)9ENT{}FqWX!6icC1VcjD@yH;IVy;!CAl8{jpojx&alTj{~3<6a3gvV{Ucmst@wQTL7EZ!2V3>1NYJ3F;m;@TLUX9bHI6(J`^|VoiiC+c^3n z+MVap0QG%s3|F_M;%jxB#Q{&LRyUPfV`7`epe9-fbQ{ptF{p+95$NxMc8o!-^arck ztH+ya|G6>T82SJ@F}9f4?lGu?Vr_1k3$$+xDpCQ^2A~6D(DHa=jL26~v}9y@=y<#v4Cs|L2d|`^UCJcN}7# z4lN5+hpMV8rM98v>!qEh%q)OXq{uJj^_91;_Xni*8oy8e8TvoO#r>RXJwINWsx_TAy2Kd;aqw9OHA9D?e{y=gIuliz68G% zIsfl$Ba}xzsdZ+hW?}nEk&nvSj#-`AReql{yT%`ok07K*Lrz!ewS~6~w;Yq&XZlk) z?hqe_VqTwo_B#w=CpwYZQ`5mHZ(l314xl1hq{xSdSkU#dS6SPeR(U_f7vAzo%b?U= za5aVFWIxl7TtM?{3`Ks!i!ton=Jy;77CPB%9y%&-#|nIJ1_OAVDTeQbW+VlZIH@wF z!?F*&{ovI)hx=-sjU%9M^|A?Gbgne{@%oEAOYLr0e)hUHY-znTU2}PhOKN+4d0yYD z)uEcu>d>8`HAm%MNIBCVfdy@e^S&sH9KTbFoRQkjQb(L+SxGYUeC)ey)a?aXFs>w? z!A8E+fu?OPL(Kvsw_IyZly;gAf+s!N*KAV;o$0O2>`|RU{%?%gKJ&yJ8j2e)JT-@G z@u@ja0pwbyR7-6In^{$goNJDeb{1?wLyml_DDrxA`lfl~99k0Z4_8Z(>xd-Zs8fu6 zWIR^RDViCThsUXn*f1JWWO&r+D)L(ycMIbl*SVXAN?UQfUOeMza9F5it zz~T}4FW)E%6Ssi*gBe`MBFLhu8Wv-(u_e(bs7_&+u}AG%U!RTPs#7@4xc{SZQ=dU& z)LB#Gx89foOsqJHXUuF4W6i#g33xjnn=gZH#Yp(8A8g^QygoLuBT{>_Uz8%BA{tgA zM$}Em82#9phn?9bKe{uoMd!}!BEKkGVlf5msjTmV(+pSUAC=$#T2X3qjqYo;gya;R z@QbAvBcSKv|0JCA(uxGXaL7yHghpYpmxdGkwy&7;R+w|V^lPs-!Suiub%n5GE6jF8 zuiWw#l_lQ8U!#4Aw+SH+T~B;Qc)?A-N}6cBs;P$QS#L58CdEXrbwlC5uF}P%jr>(w z=iVyZokV|kPqbapB0JJj0b=TPk+!F# z3Hx2NH>D!|7c6x;UFA${`C<8XWD^LfO)+(&NUqcc(K2M@izBoub&3l{&OfFKLB}eW ze%bTO)FT3afmZr{!+)_mBO>&~(9VoS!h#|CdqxdEv^%4KS|-#dTtNIDxP&BokiC$@ z<&@89bV6BDCm3?W=Qui829_U2O^*GXDkf&yEn3VcerKoFiQ9!rJ4N|PIB26T|F)PW z9J5>EFxZ)mHu7aAI(~zrD(^%&3*ni+v(kdhmAp*PW~SN?<8X&u>;Wsik@+X#WHdcA zsW)a3P))x)nhGZe;_Gm8C~r4cf~aGZBcITtlSlY-v?Xh`P;aKcW-X68hX8KPsuT{G z=+W%u{8{=iJ00)uvghJGC+8+%|Ho9H6P$E!Iy*MYBL5qzE{j|aRfQtoR!i+S*8r|8 z&u010Dm%~gyqn_^b_Y!1jrCibTQ+Z4vw6e)Ee*}L`d0b4%-mVEH~F~S++5#kZsUe^ zjko$TYvB;?mOy5v@9z5=YwOoOYo0uE`5$%% z-+7yrOY=7s*x!3$BqMnAauzt_g0GaG{&~VHNx}DCXlek<9AxrYHCJ)7V8t3#No$2h zxsmzRYS};k={+d*6}o}6;Qm4X_{m;-@Zbmc&--tWJ*e{cMCp2doet`D2WfHOoiMuu z`#6rP`%qDWsMmd@D0QG)L0egSP*GT*c+MzF7O3@OMd<>436vexThA)WNziMc>{5}1 zg{c8u4%!5|1N0HlYoJ|>|EZ$<8)z3OZq{M$q^u}s0eC)B6n4SwMcny8J;RDp1zHE% z40;lD2g64cC&F;pWAI}Utz7H@oo|a9U5Wa@0%42*Kc7<)mGXi{s z0eZpBg4>0s?zE!(nqkEb@iC#)5fCB{G32l^{_pX~;P?N4|3{VY0robYF1Vf9=Ojl^ zsI`Y3*3dX^3EDEA2CUgj@I6=J5TAr`&l6FF4qv;e(2=#v9C8F6vxFS^ZP7)J4n8W( z+a37}9D$G{3o?ZcPoYCBRINcDg-vQMTUT_Lp&}Z#1la4qvQ;dsf$ISd-T*ahE3j{Y z-K8tD8SDbK9D$MviW&HXx9kpU={U7BHt+@n&AnsfUjx<;EMJok2p>VdY#hgIu^<~l z4tQvKl7$7N_$fU-J&vb!mXOviC5OeY{y*gu-MMAY< zgf7g;wQNIm>nKf_`LmoMCV*h{ZY!AQ%Fl%w3eIJ<?;c*F ziusT8V$bpU={zic8FYh@5lWnHUEAj0#&+ ze9>|NwpquNDX24?d2&H+p7BX*NIamk4Tl?e8!9}>s@9G+3+_8R{k z4AH-;6((9y@+M-P(?U(i=z4z?JaSep9j}R`-FIYA001ohYzCd!brG)Ce1b& zEAQx5I&?P}3E`B1Pt@dYLYtNII{q9B7!<>+U+K80k625QNDx>*^h%FQt~|F|7r8z? z_SCZgHUb8Ivuctv5>=44OKv+m_m~GelIhjbN#PB8lhcL&L(}0}4I`Nv(4*tO(g$Yz zg;UQo@Sf1+2lX&82E0ed2Xu!R+3-c+B9}~eNG>lsYT1x{`i@o49oXE@2JC$oT;}TK z`mZ)Ut#Bt)?rD8I0lnX5eViXNE-v<`C3h;1)nWDr*=|KKjs=z>wMEULy=Cr&b9FhN zK32nVB|1LfQ`Na*jc`kKye~@wj30(~0iU1-%R|hr{Rv$z=U&O< zui<+N7d7(KCzyJ=7q8{<@#LzQ>gqFVqcUg;RRDb-H(yo^nnRrxQ|(5!W!br6md!}A zJIJrbE&=IProEm4Ws6i7`Xm@3<-8fi~ zivL_$;_Ag_A6rDqs7g;YQkp@HU#hBP=e~GqsuEpB{x;|UbyiKa8Gjb@>8+~S)K`W7 wpVgwx$eTsGSlxJ?_*^K=}o6N|J~G8~xxe^9-+b^rhX delta 4622 zcmZ`-3vg7`89w)Jvf0gRci9bh^Vsa>9r7Rv7?4U}Hw3Z?h)Y6fKqzP>LU{#@h?X%y z)UskK;inFE(1x}Q*v{Bsbj0yNqM)FtAho4Vm7on=1R;bVujuV}?t=;G^q;x+JO6k7 z|D1FGbIyOxZpYt*j#q>&Zb!=qrixz!*4Q>`w{@LW)7G{9t0+1vSyQr-X@oF3&}i0F z1DRdZ1pu$|mub5H?%VOGsOcwYgCIFVx>ZQG3en5LXX)+|&1N_bYBU>U^m}n;+W+gK zOYfqKtfqyovPHiw->@(8QFqR+We=ZZ$KJnHZy9 zbJ_#}CEf{-jvvp%NlOvbD0Wike$)0mqGe|1bUm!<5$p#FM0sljxN=X;%XG`DDUr7# zvKZbp#0auQQ`sW^ZIGZIR)hCB>D8$6M15LS6Zsw$)H#c$IEz%GPozp;8&a#csq#y7 zJ_?#ZD=ZQ%x5UP+8^T7=i(q$v-86(X(VdH55*!f@&4PhPB3ierREcKZzon2PHnef_bPq&Pj%visYHJ0-1-`#=2SCbgC zc`F#ZY)icXy=<6XI!rI5T*uqZb^CcJgB#D=n&?UvdneLdmnW}xB8(o+uM2K0dyf}` z;ohd$;O4RqpvZ@RR22CEboxQ$I}@qhm6_6Sv&S7sL3nNUNRDT7qIqmv;K`_CKPAx>@yJ@ZS7cf^_Y6{v{O}+xbw4QF zq+{-Cu{()8i95uMBzh|`!+Jw=RT8=D)}+{|ozs;eIC+>BL0x9>1 zGbGxYl4-rFwV5NK&AccRjO7bgsV^m4xJr_z3h(_9kO`tQLRNFscDnIdq*6imAk*eA~rA{jGO^I3p8~O5O`n50Lfkw_gts23{B3BM{zTtaM6fV(Y8M}l&@@0m_ zu{Qc^=1j5nBHhfKFI)syM!(Bi>RK7c*qZ}b{_pZvFnm%jxIiPat42CdAuqq6DC*vj zpFl~DyFjb6v!XwZRfCl9cpSB7ZxD~iQbrDmlVWLK&idFsq=&U5Fx=7$G33uph|dC1 z^7)ML8SJ(E82V-IJmDO@o9l@VWA(uG}Kx;@K#=l3#C<6@}XKYsCpsw5yJ~(#Srmp7+JAg}0DT*6d|Ba$-06nm->=)R^X+kl6F`+%FUf3E?r0sWC$2waRl8k3qD z#5WdkrXor`urg7+1JBW~6s3=md2D*~_OJ0d+(Z`_unEr=@I5-mf!jQb=Mb7iB2mSM zcRQX6g!U02n;LI_TC9rqiHqayQ{vSFmHQpQ(*eIUZUf=m2|P!7lxAfVXCRzBkB868 zP>$aWA}~kyb0Y|H6LW!E1TG86g6QvXbiAc9Ms0)7bra@n88?P-`wDPJ!QH9Z7mFuh zUloJ8AQ}QWfJq%*>4W2X)bR=4&K>b}8^!7I#ip@Bbi6$fZ-T6JsH-?{`1I*|7oD2)pnw+2ozx{p(4xtWohK%{1$21HFO_uZ z?s}Uln3$Cr&EC#Ovn7-kNHE2)T$&ZA6zsGYnvMB@QWDdD1KjT8VD5bSCNNrXlcQWp zbdMQYt(8qIf%!2)39$&2C*Vj7(Q0TR3YR}Orej#MeAS9r3q8Yl9$e1jl3G1wbcu0y zHM--wkXIYCWAGbM5Yz#@Jf%`xK92gQjK(<9Dx^f?SA}8vg<9ZaomC-K8^11$+M8-a zZf4A$!C%3Yhzw^FU%o@cfz>(!q&=Trt(Yc;lITW7wy=!SrWRFg*6WJ(?cSmLy`=M7 z^re-i^RIxnV_n%pLpnh5b=o`Cdu!bpbqA@tLaI8Yx7EP|3Q;#Q@SM&UDV+WaxH4H zkVL2O6bY45^H%5%zNU&bYwR^ucIvG(@*@mlWQSJ_-lx)xX9X_r=zL4ArW>*RMCbeT zSQ?r2gvvA8JQt z)cI9<<1P9^R*xgpRq6I!!v;26=R5QdH&Wn2@OGAz_dm10On2ZxF2?!sF?empdb#!8SB?Jk0UFT~Fm>F& z!=w4DX8?+EI^Uk9@jUWeJmC2+Vy7^HglZ|hO1Ep&17hUQc{=ZNb*5K8Ekv$dO7KauILJzggDst?J z)cQ40qGMo<)M#KXIckccjl?LMsk&yIb%~yhmr%H7{H+vir0NhIt7*nfx2@JojkR{G z-Jx}4r*Q3fYN@qj0LN;JMi>)hU_v5|m_0h$I4_Vut7lJhxUm6;iiLL0#^xUOXU8BN z3wfz;w%uW*S_A9JQCCDQb;Z W^YK=>O5;}2ENF9PYFZLP=KU|M;Um)k diff --git a/examples/08_nuclear_war/example.c b/examples/08_nuclear_war/example.c index 785a37d..d68233f 100644 --- a/examples/08_nuclear_war/example.c +++ b/examples/08_nuclear_war/example.c @@ -1,13 +1,15 @@ #include "../../squiggle.h" #include +#include #include #include double probability_of_dying_nuno(uint64_t* seed) { double first_year_russian_nuclear_weapons = 1953; - double current_year = 2023; - double laplace_probability_nuclear_exchange_next_year = sample_beta(current_year - first_year_russian_nuclear_weapons, 0, seed); + double current_year = 2022; + double laplace_probability_nuclear_exchange_year = sample_beta(1, current_year - first_year_russian_nuclear_weapons + 1, seed); + double laplace_probability_nuclear_exchange_month = 1 - pow(1-laplace_probability_nuclear_exchange_year,(1.0/12.0)) ; double london_hit_conditional_on_russia_nuclear_weapon_usage = sample_beta(7.67, 69.65, seed); // I.e., a beta distribution with a range of 0.05 to 0.16 into: https://nunosempere.com/blog/2023/03/15/fit-beta/ @@ -17,7 +19,7 @@ double probability_of_dying_nuno(uint64_t* seed) // 0.2 to 0.8, i.e., 20% to 80%, again using the previous tool double proportion_which_die_if_bomb_drops_in_london = sample_beta(10.00, 2.45, seed); // 60% to 95% - double probability_of_dying = laplace_probability_nuclear_exchange_next_year * london_hit_conditional_on_russia_nuclear_weapon_usage * informed_actor_not_able_to_escape * proportion_which_die_if_bomb_drops_in_london; + double probability_of_dying = laplace_probability_nuclear_exchange_month * london_hit_conditional_on_russia_nuclear_weapon_usage * informed_actor_not_able_to_escape * proportion_which_die_if_bomb_drops_in_london; return probability_of_dying; } @@ -32,11 +34,36 @@ double probability_of_dying_eli(uint64_t* seed) return probability_of_dying; } +double mixture(uint64_t* seed){ + double (*samplers[])(uint64_t*) = {probability_of_dying_nuno, probability_of_dying_eli}; + double weights[] = {0.5, 0.5}; + return sample_mixture(samplers, weights, 2, seed); +} + int main() { // set randomness seed uint64_t* seed = malloc(sizeof(uint64_t)); *seed = 1000; // xorshift can't start with 0 + int n = 1000 * 1000; + + + double* mixture_result = malloc(sizeof(double) * n); + for(int i=0; i