From 607554f22b67531c8b0e75cf3611257b9bfbc18f Mon Sep 17 00:00:00 2001 From: NunoSempere Date: Sun, 16 Jul 2023 14:38:12 +0200 Subject: [PATCH] add beta distribution samples!! --- scratchpad/scratchpad | Bin 21736 -> 21816 bytes scratchpad/scratchpad.c | 68 ++++++++++++++++++++++++++++++---------- 2 files changed, 51 insertions(+), 17 deletions(-) diff --git a/scratchpad/scratchpad b/scratchpad/scratchpad index 1b34cb8df77aebf313d4edc0f8b127250b1cbb85..462f9cdd009246e5d78607b0d738dd492b59d6b5 100755 GIT binary patch delta 5564 zcmaJ_3wTr29lyB{D23!C&@_G0(4`X4HVshN9Cb}->nO2M`i6fEWOb$~*FIW%ST zVZ>_j_{}{WbDP`nv3<;{V@_w%mRj~Kg4)594w2pPaG^oRDtPz%pPTS--?#gH_x{iM zzs~Rc&;R_-Nz~^I`=2*-rCGbuCGnK#MB27{#_?pEUGLfKe6pb+xn!K)Eg9Fn&`@BQ zHaJV~e~~^L_wg+Mp8VGo_THv9Ur9BY*Bx2f`-Z>n%)|@jrQO@doE(#}F+*3Z{JKt) z;&z1@-4Bd!N)9NU#T-2EVsG{0Maqt|QOIM=V$5|&QiHN=zTQbG{9=Vj87wT^~6}+U(_;9 z(`bNj4h%c(#d9&bS_QAlF4Jl7$oZvleBb1>x$yTBIDxMKh;SoZL_gPEG(lZm1bdi1 z9XEZ1OtS-TG+@l%sbk36wpe$P3nJ{l3;Q?6Lc!mep9iM86%6R*qb(3)#QLd{R1xgM zjMqV^8u6(J_BVC~hcMGnc-Z9W4Gx7>t02QhBHmz@Gp4Yap%@F)8vEct{Sl(11j9-= z;=sIG?=eyg;r(z{d?eKb$2-p%n)j+dy{u{E8{)qt#QqskF%tlN5Io_1KgBf0mKz_y z2$)YuJqSy+0Us>ppD~AuV4z(Q(TRn$CP8f+K zrh(7vrf)<}$tx*+h49#@VpN>E`6L>~|hq(pUUaL$1f0V zTJ%?hT`k-R;JL5{8!G4?gpUCl1wO+5;KzFI=##}rr!I*=v9^39kh%>$(oHPao=mJ} zJXk;(m>IefrX>liKZA{0f)B#kM<9h_RZ}H2!CZX?*v9b-$6nyKE#jdb=L$c6E4@!_w zvHCrH=oN6i8lgIYT^As=nqxrm2I?>C2q0%!P783L9&L9zl$~djJfYWMDXpvi?EuK{ z>ErJ-txY3hJO!ez{SW18(XY@Jd5&Uq3apyC$a7+t7MKux6BT*pIs%+c&NMbISyuj@2>?AV znnZ@K(IWwI&qjIzVAGnIkPAT8FOX#kJOZDHX|>a$4xxHOuV^3iK*18|fD*k1nhcO; zHQy5ta46!tI-meL>$?8ANC!4F%LxeNIEc=Z(ie+KzxtDl@x=TJvMll_rCRVkH49u3 zXdzdZp^Nd269l#^MQe-BzFM@xRP-W?kn8H`ziGv6(W2C=Cy9mw9xGRP#d$eePTiE_krd=imj7j!jS?p#T_(1mUV8cL_pN8~FPH`|n5#?On; zZVPRmpbf_}toKT$6{idRqRDgld|RgcP7$xQWzT#E7E<$#iUu|TwyuLnExI5Q0Qw_@+oG%0^^igZDWqKB z3)60}Tu33mHOB;;Nn3aFA7|v6_e>TWIdcm1_Qd$Uv=aI07#~T?O)eGc7E?(IpPK%V z+&6{qNq<!*u%iM#k<9lG0GFODTE>a?lK$`tSNMl}Ba~XnR za8JVZ02Vj@MF?>53y(n93+6wO9yrOYf$h#OBZvnYnR|$}{{}9-xHRBB5M%X!V6H=G z{z<(9Stv%A;ifIryq1c9=WXoj0kG@Cjiv0NkP;qtoVg=$@+zJ6&Ppn-D4|_u9A8vn z*Br%12_6TMUQuiPflkv?*k+>Uic?E{@?%ZATR*^X1FBJWJ?f-5;NY1gZosVV7P?FB z=%a8g?QKj9NV5YNoDzyyeLphA?<~Q$-s!-XN*KgVD0RCr1WgBRivKy0Q%Q(h-8pzU zpsURXx5-(ku7WZm2y>x^)N}Qzl{!=RT6&R*|lGPqSAM5Et!sE^WBz4#&bs_YK+ow)T z_;X-BZ`0h2bqKn$tLSkMPl%Gk^fH?IbeD&*tWmrd1JHe__%I}-QuKW^$eMx(VtbxU z{A~RP)N67gCUt`%x>Px;oqo6oHC}Vus#S3v?xze?nd^NotHvo*cMR#_{v!q?eiUaA z`1A`T#5$hG@k|*vOezDg{fN2FZXxGI`;ivVZG`kgT)B{swHL3QGX&Ldky(^=ybdor zLx$pmpse)+e1koh=`p+~4~>qnnSe3#1Y&p$@89+sns0|!%J9Ab-1Z@aA7W$Fe~3u` z9Si{uy8=s!JGgP&M#l%v*wCVc9-#XnCMf8 zyZ-a%7CIwK1zRHspYi7$*?~`RbRx0Azp{fSCTwWv)R42=%2$oFK*&9CWBGTgC)>8@p)hVwxr= zF_E$Jaj0VFn`q}NJNuHb_d~sORL4#yj_PmHFp2GW%>Y!wdbZl(&G-_b+>iCa)ABVM zY`o`f9A^Z`VtO6Y$l+0KHewSYH#q;Vx=?xX|Nah zXIv(r>enLOcTDq6c%tw`RQK2%2>H!n`6_8Z)|K@b^qNBE7q@Gb^nIuO?R(cwRr`9*K^}|X}~8*Py3#jwFo~I z@$kx7a{|wzfn!9fd|%UsKnp+6GzlwG0!mfNa$M8s*Kz2SriCEi_fJhL1b+~;6jV}i zX(WBnji6BkZ!c)!IV=-sCFmv47ElZM0d<0QVSks%(p1S2)3k@dXaRke6hJMghlMyG zYq)jx^z@4#X_^}&Nu}datR=~bcMjf-yb4kmfD1N>@ZbQHHsf{h-Lvh!ZGes6XyPHQ zXYjUyccb9w+Ak9yJn?UcVz#U@Ckg+50)OBXKRtV{ToU2Ab8fk%1o?m;qT(TMG;u$8 z0v-CqYXWZ#c(a9;B;kTBwSf2J8NO#uMnE~MX*?nI^&jog-?gRF4I#d}G)I=E^Os7O80PT_Wveqn@r;-c z%=FxdcutnB67W1yRvJ^x-B3lgh6a9Lxjp^9xQCXwheX=$;6H}6D=v-Z=^n(|qxn)2%NJmukL%wvP1hHx zjOKsh{9)V)=|QTA=X1fA_*1u9)6=l4tq4+D{BK9v1@vU{^O(X-kc#j+1XD`vpE*KW zPPTIT|KPckPgCr?QWRewVyT{$OQsyCfZ1;L>=$|v0Ewr-0< iPSAnwiaRxfVI!}oo;P8iqWk#;4_41J+{MGy)Bg+GO)}R2 delta 3966 zcmZWsdvp}l8NWNp!sfMiH?m}tCD{#KWCN09W1^2M^8`Dv^_w5>=J^cq6VrIL1m37!$X8!(U=0$-@P-@O6Q!J z@80k6yTAM0Z|>aL@v=#M$u!`!k4_`v3h09Po-8_+{=8l9J#S~(rsDK1X1%+`Jh;eI zY%1(2)B9J-%?8G1wpDs>nyI>Vy7a<{jNac>q&$1uyYJopK+TfD?vÝ|VjW8JN~ zv>0(6j=_hgwx!ESYBPW1CPeA4d^?2@d6`_lY^h4760#vmH|&%MNl=@^;(|dh1l(4L9q=Pb+VqZ^Zh%$;LH#ozf!1!IjUjjvjXjoa4U4dU1ga>JUst(r zzoT6Ty%7R<@2am{*0eWMfnM7`YXy|*02J`6Mq7{pgZm2*8X5LS)m|9Y0z7IWd#RfQ zBPvKoWAUt@UzH=na_md{Wb1RY9QF^xBN(}#Bo#U*nxu_p1}OC`X2WMc1$=u?o7V1B z_g~U9jx;%9Vg?8mofFhW><*rnByBOT1OSLf)$b6ZZpH)4bOB^*B74-e=u&r(7{tN` zv&in)DZ76ZV$AFx6+A!Bk~z6kFv$K7cw@|rMr1c1Z8v9n3YdaAAM940-_ZmNb}&dl zNsAA!yNAH)al<&WMiH_La+EklB_KUy34LQG#Rgc-PW z!!{)PgRo4&8vJ<0VO%kx_fffY9=MngvK7vL@4R0L?@d4Oxb90ve zF_YIoR*4XmDgrgf*nC319zYV)ZsbO@#KjHcoerhX{03~whG=J+7N_nWtN;O7A>sca zH}+Q?oX;-ix}>M(v)sJ=IZr`J*3Mm~^j1E|n@Z&K&P=7RvI`yczW-n?KE(`BX)(ICI%$t{U+=b9@Cp0ZV8_dVh;frw0#lFp( z<=KZWTpbSyX%BiW=;eomlg<2`Yg4vkuE7}5Ssi+`AsfuCdU%`H!ZDj%tlRYy$>d_s zx)#qnoX4TYs z;5Hl;bTZfVvK8)`)|MRp=9-T}s4e5HBR?QDkF!_vXQuCQiqvF9o$QnRN2H*WwG=!q z?aE%C5>{ZU#nK{pm)U*be0U<9e}2rW58H>N>}JRGOK??MY96G1HUG<}pY>N@~O z6`OEYkx(^^!uM!0NSopznfmr5WZHoAg|T?3ho&P~`nN^a0>~-0*^+qI8tOX)*LSGz zEo8mf>r2YC`Ck|YF&-d^Rf1^q+SGY`%LM_nxxo5uz&hG0mf+jPi6H`9fI&BGGWK)o zyBwn7_R%D5p1{c6-daIS2UfUf^9J?Nek#WIs9TbvybgkJy*G8T6nvNms;^Z21=!jMvcZb?V$K$YkSwX!h|~poo#0F#Dz`dR!Y*1K%TJgs%?zMQ*MKgRs(9 zu?mteBu`w{3!IirTlVmto&cNqK`42@a+c9(_A1nO78!*0k&0dpP)pknVar82kqZm; z?FC(j%q`Tt9wo&qs4`0yP%O{`av1s;Cr>wKq5`S! z6K)6LmXE0Kvo5}Q%zGxA8;1FM9Rcbd+~dZ44`aSVFx2lw4UN)`yo!6HrpgfzpzG(O z7uk);j^I9eU?LHwbN-p6b9niU@RHpBBKGv>1z9XsDr#d0;6niY5awq1HNK59Y=oAl zw$q5kD8Dn528 zrPJ*{YV7obqgR)v4c4V~m^YP|KGrmzz@ASu?HuF)qy;l-J*H_7;r*@~)3oi-2R_p@p43*za>xNl z1#%R!6|(0GP1^z)RWWd^B6Wm#ahKyBWv{^cIEV^0H{rr|L^^Tzu=CgER5?Kt2M%#@uWh)!(9IPZ z-TSX-2cQ!NYQ<{1$7&Jy`=MVw#uk-dC&iDmUzOi{a~zYwxkp?atpO7`1D#-pKG(HD zXE~v1*9tE}@Btz{(3PKHm&^0QDoj26q`A;)>y;K;w@WECDb^iQtu-pGfLsJkOS;u2 zTP<>i6frETA!S$yxe1}i0T=&B(@vpJs;#!)NK35F{*)@K=hvyS)%TcLwwCrzNeM~m zR^I}vr^@QAvf8Swmg)?WL2xc4UIfNZvr`pSmNqmkXEcqS4_q}^I{T%m_Q@ZE4t!Pk zWnwGvdF-h4NLvaAH&x!23dzRgVp9*ZSLK+Ttf=bWQr&bWEhsJ{MvNXw2v359w3J?< zuQv1@4z_YZwY1&AUWV9Q!4eC)fO}tXw&@S-$)H;rC}gh&mz#=N+QNJD2)_ONx{1#l z*IeWwTe1`Id17I))aqip78ct_ah&1TZe5T|EEdekoq?EaFo_#XCJW#iCRgW}Z-CXq zUXiXcDa>A-Q*#&kb}5%w>I@J=tTs@h2Fhe<@i)Q_d{U}A@&)GCMvj*Iv=QYg)V25M zg1pSFtn?OVxXkdF{1Kcis5cCKz?joy@r^)lBYDJj1-}mr|7|!?@=H`28Vp24kmXeLvDZLPmPtw+%()yF=pzuuI=?8|2V7pU** Aq5uE@ diff --git a/scratchpad/scratchpad.c b/scratchpad/scratchpad.c index 764c783..7b2f4ae 100644 --- a/scratchpad/scratchpad.c +++ b/scratchpad/scratchpad.c @@ -211,11 +211,11 @@ struct box incbeta(float a, float b, float x) { if (x < 0.0 || x > 1.0){ if(EXIT_ON_ERROR){ - printf("x out of bounds, in function incbeta, in %s (%d)", __FILE__, __LINE__); + printf("x = %f, x out of bounds [0, 1], in function incbeta, in %s (%d)", __FILE__, __LINE__); exit(1); }else{ char error_msg[200]; - snprintf(error_msg, 200, "x out of bounds, in function incbeta, in %s (%d)", __FILE__, __LINE__); + snprintf(error_msg, 200, "x = %f, x out of bounds [0, 1], in function incbeta, in %s (%d)", x, __FILE__, __LINE__); result.empty = 1; result.error_msg = error_msg; return result; @@ -286,8 +286,16 @@ struct box incbeta(float a, float b, float x) { } struct box cdf_beta(float x){ - float successes = 1, failures = (2023-1945); - return incbeta(successes, failures, x); + if(x < 0){ + struct box result = { .empty = 0, .content = 0}; + return result; + } else if(x > 1){ + struct box result = { .empty = 0, .content = 1}; + return result; + } else { + float successes = 1, failures = (2023-1945); + return incbeta(successes, failures, x); + } } float cdf_dangerous_beta(float x){ @@ -295,20 +303,32 @@ float cdf_dangerous_beta(float x){ // But it will propagate through the code // So it doesn't feel like a great architectural choice; // I prefer my choice of setting a variable which will determine whether to exit on failure or not. - float successes = 1, failures = (2023-1945); - struct box result = incbeta(successes, failures, x); - if(result.empty){ - printf("%s", result.error_msg); - exit(1); - }else{ - return result.content; + // Ok, so the proper thing to do would be to refactor inverse_cdf + // but, I could also use a GOTO? + // Ok, alternatives are: + // - Refactor inverse_cdf to take a box, take the small complexity + penalty. Add a helper + // - Duplicate the code, have a refactored inverse_cdf as well as a normal cdf + // - Do something hacky + // a. dangerous beta, which exits + // b. clever & hacky go-to statements + // i. They actually look fun to implement + // ii. But they would be hard for others to use. + if(x < 0){ + return 0; + } else if(x > 1){ + return 1; + } else { + float successes = 100, failures = 100; + struct box result = incbeta(successes, failures, x); + if(result.empty){ + printf("%s\n", result.error_msg); + exit(1); + return 1; + }else{ + return result.content; + } } } -struct box dangerous_beta_sampler(uint32_t* seed) - // Think through what to do to feed the incbeta box into -{ - return sampler(cdf_dangerous_beta, seed); -} int main() { @@ -373,6 +393,20 @@ int main() clock_t end_2 = clock(); float time_spent_2 = (float)(end_2 - begin_2) / CLOCKS_PER_SEC; printf("Time spent: %f", time_spent_2); - + + // Get some beta samples + clock_t begin_3 = clock(); + printf("\n\nGetting some samples from box sampler_dangerous_beta\n"); + for (int i = 0; i < n; i++) { + struct box sample = sampler(cdf_dangerous_beta, seed); + if (sample.empty) { + printf("Error in sampler function"); + } else { + printf("%f\n", sample.content); + } + } + clock_t end_3 = clock(); + float time_spent_3 = (float)(end_3 - begin_3) / CLOCKS_PER_SEC; + printf("Time spent: %f", time_spent_3); return 0; }