From fc17561028ba7f1106e807788484077e1c4763e2 Mon Sep 17 00:00:00 2001 From: NunoSempere Date: Fri, 11 Aug 2023 14:26:29 +0200 Subject: [PATCH] continue nuclear recovery modelling. --- examples/08_nuclear_war/example | Bin 22448 -> 22672 bytes examples/10_nuclear_recovery/example | Bin 0 -> 22944 bytes examples/10_nuclear_recovery/example.c | 63 +++++++++++++++++++++++-- squiggle.c | 1 + 4 files changed, 61 insertions(+), 3 deletions(-) create mode 100755 examples/10_nuclear_recovery/example diff --git a/examples/08_nuclear_war/example b/examples/08_nuclear_war/example index 47874565c39358521080f3b9ddceabe77de8f5fb..8df8418cbb9c1e2bb2e6352217827f24f1636b91 100755 GIT binary patch delta 2562 zcmbVNeNa@_6@T}^$}W|C>+-RDD7%t{EFxdBAZ==gHM{V@4N8P6wSuUSP1yweia28n zlNg28O>s^#(Nxn+`Z2bXSj{#qX{|+YrljpmFiktwYTQKAY?KzVng&hu_IG#RwrSFT z+B@_1ch2ved+xdCocH$lPuTDYHtdimR~=aNX7b}!DY_oFl9`nxSL4$Ik8BW^V@4Xiaz|~`rz(_Mn%q>N540Br{`x&lKQRP_RLc? z%N40hRlYGEev7=R&lvXRz!FGNa)Q&TPs^+;xG&>KvbRMYx~)DoWr`?kb$&&XBFenE z_VALAIai0Kp}6h(h+@(;MNMlH?O_W{V7zv{`q&xMdh9zLQQ9WD6vt_r%t#B`7dFex zPP^{<0s9hJ?cZdp>Bsh!^5tvvk$s&!XNqk0W%N}3gK1ABOVTShaJy6DM;*EpK6#bQ z1@#%_Fofq`ft;{;TNycmV$i44-FziQ~I;?-GWEu>jt{(Zenfp^Si%ZVuI-r zGmMo+NL@57JXbZOlM7iq_}*n&SG<_D&|vWw`y)#0tLV~XmohyH z`7e$E`9UjHBo1Eq>GN;Ukx$y-u_LtKq_B3|p;xqEmTY z$-RYmW8;c-KoBwgKhEp_Uk)t(e-3o5@dlO^EXMw#y$2*+w0j>oe1P%gU5O(0Gu)oXs1>iy8j4bs_=g_FR2{_;? z=MLZq#@s9PY0h^8eZa%O3E)X!^$h2ez&_v=;4E-f@Da|_FbuvMoNoX+KEOVqe+Yk| zdW-WK%;GF?Gw|p~oDb^C;D8r-0>lI`1Uv`40j$1_#4y7u&;}d<`sUEsl7j3}6s1)3 z*_~t(V8`RcGkbyLih``ki=0P!VZ%zpvmMU_xhgD}&(;bn8>LtIh9;@TGWk-nQD6T#$p%$?cfx;P!*t6=l$G;2ChUQ=GO{ zrt3nG=To4vZq8zKSDK;ol~p<)#66T=RjLysF--MU-MYDGu0$8Batt;IWzhCRv+9~5 zWfEJuv|rT!NoWWE!TB1yK31?uyFZrnCz^k$i2_vvvt62m3vF+CHET;=V7W#dmG?MKCEg3?_w^Q3sVL)7Ur}4s(}g9Q`L0<3bVQmSF!qD6)5u=5BFYDFB|lG|h^`5h zI^z|v{JmVA0`|=ZItZGA-KE>*gzQX1XXOWIZBZGYs17v0V$9Z;C@k++o|X|(KJ4zi z%lj+8f(r#}Xc*=g=E~1we=tnh4>g89Isz{Z{T>f$vBlhI;bIf211E7Y-oGw)JHIa~ z^LcPlmCj1O7xp+^23AlWcKemSL_*wtX&qv^vPdULioRR45GUT7!x&TPYJ~m77kXRJ zrElu%1^prBj6`6vXan5}-r@#H*X5Gh$)WLxd_TN-tQ*xqq<{$Zh4D?sgp<>bGFE_;-DWI^JcH~*w%lQC z1UZG|IYJ6>>+6sW$nD63$o<$Xr`;3tULi+F`D89Ry$B6|e=F=H2zvZst`v0p<-FH-h5&KLUkoUtin@E^vqrFbU) zVxdUc3oKBIx1|M*At~HnkWV9!TA|_e*PyT>pH9Kb3QgwtHO9uVKk+9UjV^qB&{3ho zfD`x)!ug6U+631s>gX9LH>F#C!fqMF5H<3B?9k2sX)UO|Ky8yaPmAnAtpPerSs6{( z&&Jn01LzroUQ-s0!mp+*PJI{dn{sFuC@YP004$Z7jY1+Hb~1_`+$@F)S}SUNlrht5 zS~qIVsC^`ALh@6n-9&8>rXmS*Xo~*b6zmV-@eQL-0>4)}6O*ymkuk`rs*>_(Ho^KT zqqG&R76?=|8lr+5kDv#Gh$BVa#W_>TxLt}AwLL}VoI&&`tMNFRfts_xYL-Y%I3nOO zn-w)^G~=lK8H}2j$Oh0n^NcZgY)!`tdsqt97M`mv4ehmTkRV~y=Ac6%!;)Vp_qiD+ zmL8)EAyA{^%1YsE%^+O?JM8ZX^M>s^@%M2#cA<3Mu`9EC9wx#k-jJ(UTpCf6ZdmEm zaL0>an{z{IxFj|y-bpJ2oEqu{Id9C_j716oM4WRjHkGuCdN>c(@cHx*bnu$;!=iXu z%o=a`Db&?uHc_t-qeb+JRU@0{l_#Jdh3m8LknVJ moR=k`iTYQR8<)ex>MBLqvY7j%1Qs+@$%j3$;cVz#rTY)C5-lYF diff --git a/examples/10_nuclear_recovery/example b/examples/10_nuclear_recovery/example new file mode 100755 index 0000000000000000000000000000000000000000..d33589aa250be81f49ff661e85990ef03eb21d26 GIT binary patch literal 22944 zcmeHP4S1B*mHzTG5FnXvKq&DC9d&5q50eBCiOM>I34F0p0){rh$S@=`kUIHsG6O+f zj17My#>QQ{+W*w{+1lOhQ`&8xqV4J=fCQy3SZ%fH5&`QsAQ0mp6=e23_vf2$zGNzO z`|R`Vv)m_{bMLw5oO|xM=bn4NduP7;DykM_XJsieiIhWyhy@vzT^7Kk}I?NOJ96;hLq|7MrB92xhlo}EwzD#M`I9?;6;^1ZC zvr$%*t+3(}PA8$PH`*)V^@;^niw;R><|Et8oyYm+-jWWbT?U1HBp(%tPZNBa;FD01 z^a4tIl1=D;uCS*SI#EZ-$vzdu%2zA+tgu|rOGs@{y6vq<<7a8FTJUWv<#MdBiVTh< zp=@t4_{c72o8&Mrk$PqrhXG+n2}Sctaj0qKoY}>p#u=femT1R}j=6JY%$ebz*ePbNtj)F#nT zB6>N*RHi)1AsXKQ`;p@~g8Kt6G!x8}evFDqqf{QqfS&-oa}@dqGth6(fM1;f|8NF; zO$K~J2K;X`;J*O86aUs!3&d!4csv9C?hN=)2K?p>_&LBk@o!Dbff!An>J0d&fxpZu zgv#*KW^&6pim$Tf7GGl^99Y%V9tnhNZn-Jc+7hVouM7ngpKn!jYm2Wv;txlBK2AC+ zVW6W;3AL_L+V2WSg6xxuL8YyAt^> z8tx>lFX(RyDb4;+sI`IcEhdJ$+FQdBAR9XTzF<>}Kh$(LX-VRRmd$=N7Yqji%EGG3 z`8WAyx@WqtN#(CiWy_50Y$G?zJsZZ!X7ecrpShT+xtKvTzsBM-n|)F~m*KMTnT0x< z(;oasFKEI_GFCY(aGM|gCrrOF%8)3>SUFcYg)$}LZv!NLN|T$K3Me{^fl$J6c{zr> z$oV8cTX9gMNOFCUak+B}$9qKlmw1oB=T-B1iI?S97I8eq3hJv&gC_iVi3HS5c$)h% z4Vmy0B@*zY2~TlJrWFD|76{QLR<+uMr&uOan+Z>GP^JzOo~_rSqSJ&ouLE05cxqdw z9ur>Li0U3Q;l);gvGkho`BDhr6DB+c)JWS*c)NjBl)M&OgNaXG zb0yHj$>~=CsFx^2znW@>6K^Jy$sSgw1%EiP-zd{UKb+_@%Cr#~PW;R$(}F*o_`Xr5 zg?>2kuu-N3emHTzQKp4`IB}OzrUiXC;Wx^(kPjyo8)e$y3@5xsnKm}ViBh9X3;b~6 zQlm@@`*32SQKkibIFTjG=>2d1j(T6Nu>PYw-BO-rDNnJK3oK=urThsVPtEc_u#}Hk z%5Pc9Z&=DNTgtz+ly_Om|6wWr(p>f~eTFRVtyy{-T%}Oo0>fzC^Y%=GzwB0YjII{j zrDqd1b{rLS_Mq>$0O;;(_JZgc@{gc;Voy+0*Y^e0_2sVo?c^iC#kIbK^F4JLP(ih2 zKXT9Mb0O&&WzyY>elc>Dv0YkhM^$$(VS3=F<1VHKPcY`L9eEJiGN4b!m<_8FZ#r33 zb}0JCs@RFAJ@ZvXFTfs*HNTy7fIfZDNB_ea)`42ef06P`?fp2}P z-w#P#?N!L7sU38ymu!FnP+_aB4zi)YCUpnuR-mXKalHz@$of%P_ImV1qaPLMURpo? zjDs;<<^9-B`L0i7KVC&w_44yscW*&en)&-jnEy!1{0GsGz50`o#QX6*U;RIW?@-*&?Vk#v!HuMszDp|1#2f5!3mE{;Si_y=!M-Y5ev(y!GBYcBJe@@ul=w zWo-X)Ep|}rIy_WeQ`K!-!$zqf$vQc&7Gjc9CS9*6@Ats{%8=2zZ8Hg!os3L@%Ezxm z9-(bqa!|dblg43(fS@18=paK#xdHl+0(zftQunWJSrtL;A%nG6BDgn+!{V35sI z$xr%q&ZBB2deu^`Yu@9O)nac&^0lq={uu@R$y3Rs-Ut<+wS5QWtUJJKv`_qrwWt*O zP$uz9u%*MMiDgnRM3#|Ke_--Q|L~oC)&fP0<;lKyo~^v)?Te` zFRQOD;Ba}Lxqm^mXkb{aFW3$mnstdES^4WRt2xahqMduGo%=wSQY97NREG*(#%xEo zlR>q(Z|%6!IN*T#XncpMuj(Zp7H_wX)nffx-@yW{v|o$u)B5y@GB;?`O83Ra#^Ppu zwUG>C4F<68(_&*Y(W4Gc9RWSh9TfB(pQD}+YwDG8vh;Lx4He<4EEbFQAJNA3 z6Gi{9TFa=&KhrzJ`a?m~tYd~YQ(fMVAPJ40(_+ZI35}4Os&3f=PeaRlsLkhK&?#R- zd(aN~iB+yVIJ~;$K~7_6-mic2&&ecp;VU3!-0HHpTJ(zQ9(Ye({Mvi!#Clim+BFbY z?t+)q*=`1%{!1{?5K~84w9Z9UAI;evd_F+6u7hNUJ=*?+Hf}%Zv8< zpCYsgE0(()7+_}3LhgVW+WyqBdc-T%2`4QhmCyTMSMDl=6{tsDJBwfou^wd{w=dJU^^hM8+ z&tlLidc%vC+cx3!5c7w5Hcjq=T{y!7?TCm63V5X{>W&tI&QrI1$LiMkXaRzG2y!sL z{>BKup2~VtXL}p1XBd}@)Xn#SkA)!zIrk3G>TWQm4bDY|bE7-8xi0v&^Y+BGAd~BR zM*qI`Q~Y}g`FHpG7@xHk|Gte>;Nw$*Wy{us13LAqo!CCoDtRRoh1xcP=-g{uLj>fa zpz3~wSclccd(_QCjNUPA3AiyfufpW@uEN9*(QE7v*zoeTq|Wv$MuKR>W4dlJsO?eI z>;3O>Xa5HjzU2y6k|r8UFE^g5r5uMsH3uf$ z4;}ZXslyN>*Ky2{i8Zb~b;~cf4cXd>-qjV{#(G!u3i32azXQZPsl&8!O!+BJz>&B4 zK!YL1+$-L{b~&h^C(MuY(^M%I1avalAR1%1}LBsZmqs_#w(aWA zM%-YM$Vktk)huf0m|10JU+7y)`ZoUc`g%TmHY+k#C=BJWCAsni*v5EXARH-e;Hgj= zr`;9n7wHcaqSKYfH-4x-^bYz8UDOT7kckbW>M(4${|IcoQ`nIC{ioQm>{)fg!;l~} zgGRfkhwDgU(GdP$7g8sGXh>!H1jYCS)%q8LY8`D^Uk$2L7$oZ*Fox+olI@thn*snF z|3y$(KzsD20yAYd$hacvoLw9!;(=!#xH`aEfpSc9pV8rQ?DW+!!(QSD4DKs`Ma zpaSuM0+il(oH;$RaBr$0shIy{Kkt{lU#g6~;9Xi7`@~!0jlG0(nb>dEzzpOnWew#A{i@pJLDQ*dQFUylhb2a3dsei|)9k;WRWu9INbmn}TSr zt16b%*B@hAo(IEV?TzgM1PS^K3r%!XZXIKU7xCbh{v77*I-Rxt@4ekMt~2a#Kx014 z#s!k{PeITC%#7$hvE~yz+3toze;hqr85@I*MaRA9S}K13BMu-+KLqvjNAqESt?St= zZ`u3nzW_nyAK*M7_XD^GtzUg>4pqfI@+NC+%F&5EkQMnk=z+Q^Su^jX3^DrU@35}1 zEkz+emhdk3`sDmH!rh2C+)C?`@w_v--4<}+pVe-OJy;<4a~BOSjcej|hC1;YbEzGP z3KV60^43(w_S~l3|4tqj|3+7l7ZFwNM&txkDr5C7N9{6gE0@@)&Eic&i?#cA zQlnL|!P~v{JJ=wX$NH*b&*3QUE(q7Fj^eP1gy}1h(~m*CDcS{DydwQiOCR@OK&M>{ zR=_F0XOPx1Q@6@pb4%m$d4y&TxeqP5a9>9f^VAU{ABCqPDuWlZt5_!+deldEC4a{j zVWXaJfE?ggH_=_aXTvdd!&qd%a1cTYx^VdX6d3`ocL;G|FXF-T>bK%*MgNAly78N+ z0cOAYt$wwfVd$y76!DLCdk{if2(~z5e^!J%S|T0m(I<>;YUyUw#V6Ixe`FGi5r?=| zM_9S&AY@?=S^Fj_bsWPYbV19|W!HvybQ1J<)zD+Z-SATyNwBk#WVh18UKW0|tuUf7 zjOxapfd+BNL*s}OI@9?QT8ogT#qX}EN9*v#;wc(TvSt!$u33N_Fz>0GBk)J>b~2pS zU9b`AmEv&lZFs>}lfYP5fiZ*B)Xg_9AYg{NmBagFWa2^&GZRADBd8 zbz?rwH_nB~%{bpoe1zuoUOky)@p7rjgYJK)gKvSSyozfDGSnu9pv zs%B(FjB2cc8_&QlDAl`a@fo5mH%?#X!g6n@oBG-K7VGM!x!mnA!nNdjwmft|xLzHW z6Zb&M?WQnN_da8|H|DS=$jP`aRcuB;S?C2{gI-lkFXSPc;h~*#BUVCY?TulRj@Zz` zcj=L#_>{U?=NcYRH+{%*Mf+$xNXx`h8X~eADD^dW(4s>14pro#@XcbkC4;fzmwcW_&uFWt3`wOA7H&(%Ytkt=+FoE0ri0(AvUhSr=TNC z5ba=s?ipxZGJ1GDdYG6d?UTZL2yYvSUrzXAkfIN~MWQxPvwS~5@ex*YYFo#f92X8< z8rLw%Vgp*=AkNMDSWd?p%yWY{NkndZ$ITc^c5a4XRA;Nib)3GFcE^YU@ZQ#ZS7i)2 zga(`m09R-jDS8#6f)m!En83yvO?15}&O>!=+<-WZtBoK1q3(d*Y(D;N_5DUJ^nD2j zhkt8*R~#1le#y%t$TP|{`3jyUKA?Lrhu(?$%Gk5)cIpYLqVeSyJKp3=m{)Wgbu{|U zQwT6{p^u*K-IPV)?QVn+6oqF`tSLZGUV672%!w*FpUr zi!k^^Jw5Z>eh2Py2I#&7n-OH=-sRZaI<;~8*i}+$zlcGl{VaTu_T_ark)S=Un|vAj z)5Uw%+K}ttfj-fKL_1Ix$#Ax$zkzlU0O*v7!ke)d#bBrW5AtYt+A00z+{$^lepF(DP)TuD%H;>H8!na!KiHLf@*`a8>Mh zdF*X(a`FLO$zg~2M)U~!nP=MbnRX=PzjQNBkYo}imH)4%EPMx%(&d*s=LdqV;ea!; zI^Ya8h1(;}mS{sL;14?+T0_^u{X`x4RgB?RQ0+R;>yJ+zr?XiT+R@Blret zQ_CvW1AK9IRUqO-Pon{6OKZd#Y>l=wGGW1J3-oGjf#5ZPaC^Yl&=^z-XO!2Rx^$e> z<9e=-66g5yJ-`0Ss&PvekM}HP>1CnkWuE3;XMTRqaIYv2Jp;bMWb!|E3%D@r!kp{r zyYWRpwEb5y*@=AHfn@Rtz)rw8mBCSiVI4Y{Os+W^Kn3yY9l#R6T+~+s76P^b zUJcj-=mCra_5cn6{t(cKxJ~!<)qs1c9&nKAvGAS-T=8}?>A;;o?YFK5986%E0d`_X zAEGjzWxowraw3`h4PX!8>wwBoGC2h30n9}Zj{_D0;_D{K*?=WDeC-4r0{jh?KLS0V zEKnI`)uoG|#;gN7i%Tm_e zOO&h*N7jY;WAl1oBf{f~TIoc{C?<;Xv226vLUuO(_1Dn`+0^b>U@yE`ov=2qQ@Qbi z>#m&Tx(ryxUk@061GXWANF+Cc|LwpJNhBcA&|RH^F%w_$Wion!$_MWARHxY;ci8h7 zEq(9a1w1b2SR($t_%E^G27$Z5f_oo04Y*4IDb2-y9`zaiv?<|Uj>aP zrzxA$ESQk}D9?Fu1yBA`4f-44_Z7$!{U+Gac_$Od>goJFg!;wi`n9~CY((?!-%uZ* zdL+Wp9L&Rs0r?#z!qI%34ct10;dZ0B*$W))ER~#TS$1a~8zT)x9d$c#T0`@k#`kO4 z^X-KXIaLl$VXS7nU<$fq_au-Ui;QG0!Q^7`&Uv z&$@CK+E?V#v%onEoU_3H8w(h}&k~z8T5>5_VH-y)`2w=Qi9Eg!$?Td4qJTIKQzZsI0s{8Is>JI^D%a0C*YufLjo$elc(emkeAN>Q~rK%;Y~MP=PX*fGTIV} zI?LR%-6bGWd^KTpSxRpfj~$LA^1Ptx(@6zN~-`09nCpmahlW`@T9#OH{jN$nlL_^1@co{WQJQSXCm3NaV49`)# z8x4=I0Kz*R%03y;1S7>I!dnC5RRS-6mn+B>FEUvsJw=#g-N_ooK5K;>xvoe#T};kM zJ3kKmX!>XB*H4(7$x0PHz**90a@aEqIxE*3$+8!jEM??4zRUZvpQNWVJ8tJ%yNLk^54EXObyhCZbn)7TEtaP41`m~YbAzjLW5(nOacCCKCo6!p*6ZjVx zgyg`rSd#mR3yB_I3A4t5KQK9n8y5Og8RSoe9Y}tUm}2bA1T9@E@D9-+J7)pDSm1lv z!Ih$TP^L7W;d$QtxkS*niH!q0`vE~f;5`=mJeEPuk2By0r2NZxtL%&f*%yH)J6p&1 z0Ps#un)Z1d_=(x)D{{S*m8S>~3(2^OM{JhH!=K4++b-qh*#bXR;xFe2b}obL=LJ4a zfdxtW$sB>-C+r~4)x1)U*suw6D_>=JF2nr!bwOY4=2&*NKU&V3_HIeJ+1SrT+lnkddkj2!Bffbg4qNbfTwY3EXFTcQhgq;I|Hp=H)u*7>yg|Fh)a(Zb?I&GukE3doNdrRd_={4*fFF@gqLcR)3NYKic zD87YNi{^W)e2W$=xUHhbSL2;uRY6kOJ69UoqfBd*&F^XPc;w4h=!GiuHkK5IUdPhT zjOknDZ*KPg#TU@fyJ6B<*&ATeQCwP>F3@26D80Cb%5`{G~Z4 zq`&9K*Vx|bTkUUYq;~*SE&^3!Q;RR!9%vL28luhZ(Pr6+=BAEFG#n5`J`gKgJCqK1 z|LUe-WX|lfv9zUczO>t0{$?RdXF{t2E5m*oKDy)e!NlB>Gzks&q?sb#CX~)DoY3G> za8v6mE5XRiYii6Keod{QiH1OXyDtz9w}y>=&6=tBg$}TW7{+#CBD!sc9xK7c=N5Z2 zkdTCNR=4>2iX35N1PzSGf5^>>dn9{w?CxW0H zcK9x?{+s{8YoX|MQ0b}~_G<8l0;v&`j+yeQl>y9G^E@JlN`LWFx@0+={q6 z=U2{%_*W^6e6_!Qwc>7E*McTEjD&d&-S*?&Kds_fV!2%Hg{mPSRqz7LYgdb8GNbW95Rol{LM`e#@XShiW_mW8B32` z^*%GC?S2;iu^nScp65s?&ml!Mqm_Arm*by~=p?^9_mR+vI&+fyCyASbjG4a{`$$R> zx_Hs}#P%mGRib3(m*+?lI@74xb~Ft``$$}IwMzZv`KpAneR;lW)!z?(dM20r^88i8 zlVn^Z(vh-~U!H?D0Y_by{PO%tLV11*5|)HGSyb*qj?P9UzdYxX@Jc}`^_O%Ku17tc zwaT(Q50kJ{@Jky@{RLX-Mux^j^2>8G3FCrLwlDc*{~s0n^97weUz1RtKNCOkSoyyX z9N9$5$0G$#0*c33fgCyh-$sp@-zoScv^u;+yz(O}|3u!DJl~V_QUO_(aGQl+p8rWW zBpC&id=frs;g`=f64Hw<%t`V|9DM=F+<$qFAYp+d6wo^Vb_sqt{_^=)LOTC2Cpmr+ z_giF$MDok?NeSh9LM6Sm{nsq~#yO@q&#cCCf;m~+e+!j%{EM(-oOjCix?1UFS;BY0 ztD5-bd75V)kVuk_mr`Yk{{$7}yRvQkCNn*5 zZRbVfQ}Rnl?=Lj-8|Ne*&LZ290wkX-=YyGWvVHly>hN$1D?PPkPIT@^xDqoGdC4!1 zNsL;_Z-&x$3?gJa9$wleF2ZOU&>4f&N5&s|ZkM`q4*ofGGl##o_>XK?($aJvjsMwd MPH};S!BSNI8_2cRjsO4v literal 0 HcmV?d00001 diff --git a/examples/10_nuclear_recovery/example.c b/examples/10_nuclear_recovery/example.c index 88a4746..670355d 100644 --- a/examples/10_nuclear_recovery/example.c +++ b/examples/10_nuclear_recovery/example.c @@ -8,13 +8,37 @@ double laplace(double successes, double trials){ return (successes + 1)/(trials + 2); } -double yearly_probability_nuclear_apocalypse(double year, uint64_t* seed) +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. +} +double yearly_probability_nuclear_collapse_2023(uint64_t* seed){ + return yearly_probability_nuclear_collapse(2023, seed); } +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); + return sample_laplace(successes, failures, seed); +} +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; +} int main() { @@ -22,7 +46,40 @@ int main() uint64_t* seed = malloc(sizeof(uint64_t)); *seed = 1000; // xorshift can't start with 0 - int n = 1000 * 1000; + int num_samples = 1000000; + + // Before a first nuclear collapse + printf("## Before the first nuclear collapse\n"); + struct c_i c_i_90_2023 = get_90_confidence_interval(yearly_probability_nuclear_collapse_2023, seed); + printf("90%% confidence interval: [%f, %f]\n", c_i_90_2023.low, c_i_90_2023.high); + + double* yearly_probability_nuclear_collapse_2023_samples = malloc(sizeof(double) * num_samples); + 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)); + + // After the first nuclear collapse + printf("\n## After the first nuclear collapse\n"); + struct c_i c_i_90_2070 = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_example, seed); + printf("90%% confidence interval: [%f, %f]\n", c_i_90_2070.low, c_i_90_2070.high); + + double* yearly_probability_nuclear_collapse_after_recovery_samples = malloc(sizeof(double) * num_samples); + 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)); + + // After the first nuclear collapse (antiinductive) + printf("\n## After the first nuclear collapse (antiinductive)\n"); + struct c_i c_i_90_antiinductive = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_antiinductive, seed); + printf("90%% confidence interval: [%f, %f]\n", c_i_90_antiinductive.low, c_i_90_antiinductive.high); + + double* yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples = malloc(sizeof(double) * num_samples); + 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(seed); } diff --git a/squiggle.c b/squiggle.c index 0150b6a..adb50e8 100644 --- a/squiggle.c +++ b/squiggle.c @@ -135,6 +135,7 @@ double sample_beta(double a, double b, uint64_t* seed) } double sample_laplace(double successes, double failures, uint64_t* seed){ + // see return sample_beta(successes + 1, failures + 1, seed); }