From e94774ae3aa05fe84470cb092d256b18fd42e464 Mon Sep 17 00:00:00 2001 From: NunoSempere Date: Sun, 16 Jul 2023 13:57:02 +0200 Subject: [PATCH] add dangerous beta sampler. --- scratchpad/scratchpad | Bin 17448 -> 21736 bytes scratchpad/scratchpad.c | 154 +++++++++++++++++++++++++++++++++++++--- 2 files changed, 143 insertions(+), 11 deletions(-) diff --git a/scratchpad/scratchpad b/scratchpad/scratchpad index 1516d8fc829a7f59f5ae97e2f8b23ee4122eadb8..1b34cb8df77aebf313d4edc0f8b127250b1cbb85 100755 GIT binary patch literal 21736 zcmeHPeRNaDl^lx~ZEz66lCU?^EBQt5z{cD$AT|*5uGZ09~mTX3_$>s;WL5TjLb-<@SBYPYg6D4 zO@VKi0zWbZ{tv)y_!&bBfXVFfO@V)N3VhoXcz6nY=M?y2;5Ph>p&7trcK!(X^+sk& zqu+H=H(ShnjV-HvfnYqiC7ehG<1MRKL}Jlki+^(@$b7ynU9qSyk@UxtKA*-qDPwSZ zH;ZiXcXjz$B({Ym?v5uzEEo^5L^w)ib%?Tj-SKcV8S-^>-o=93!%5cNlT4r(l+bY( zp}vqm9AO=iSjSzg%O3%y+^9~>=j+(+_l3eye7J{@STRg2i!=q2rvWABk6GbxIgRow0@L_QI6v_uKSJ`wtdt5xkn_8o zUmvg3@MbZOOT1g)rPM(XBrem(*JwD6c{!&LGOEL+4B=@V4uf*xf)1B;qRfjrT;D&t z1wI`-qDoxXb{(#tpL%sTwS)|x)!|}^U%ds}mtBL&Vs*sZcc65bS3s>@1~{wEy%|KA!-&6=apEpKQjTOVfrY zo4Gxgrj0^2nXHXRdQlg{Y{@ULm`V6%)zf2A(j{RfG`lD$dJV>e^ z)l?5ab^3c~Jnx?J9|I!YEzr$%dX?x9DduXRHjNT~47JeyB@nI{N?%CN5K31~S3gLw%}? zjd_(D2FH{(^B18heIB(ZFHD2oeP@e$o=P9Rtg56)7MR5iuBLexxaC~-Npka;veCR5 z3SjqO`u7k?Ka4-L-g~gMF?Ar_gCuY7JiAi+ht5>EMnF!RT_fDex5~L~XOQ7>{k%SP z23d%t_ETYK#WwR9sN|MYS*tbkPG*r&0RyrC73cA7L9+gfX^WPfvwjh`el*|k??NFr z+zsvxuE9G!?VdY58;7~QuED0%uXyFt*C88fm0kPb<;-dF_xbl#_3XoFlP6}6x5-l| z1Pzi!o+eldLGauQp0g(qX{5-pP_F;G-f2u7+vrW5_V&Gbp}D20-~26%cyE7sH79jx zXi_UrFJtUC1MulpH&dp$n8?zL$x2LxKSdr*g;%Mg%A8&_2!r*Ug8x;@?5}-^$h<4w z)CKR*`CGh07fs%xXT8U+BxgavzX%1T>ZoRyl&AT!w~n$|&yDN6eRXk)dQ)#D%e;H) zwjz-Z{{fx|z{F?GT}T(*3E5ovB%j%aOkoeD6~`g!QcCsEwlW%yC2-i7a_11qQ9cbj zrTWmea%7gjJsUGc_0ZO7ieqRvGyqw2fqt<~ zAWBOF^^W0C$f?ZPiAwPGVXgQ6a~BZBd{TZ0Y$&pFJEs-(c4)O;6XHs&=C6YWyn9OK zI>7Et?X;&sU4x=tD6YJlbG8ZEYTm^4SR-~6+-6?SbJwG0FYsRaz?U#`lV9@V0s#QVMps5sOV9}c|>`1P-)n|6T|WeWGK$(l}Gn0^&EJNOgKr-98)UV?N)eU z2(c!68TsP))ug8)VE-uoTkKWo*MH9kGxyCf&$M#ZHuF!k6jYFHV@f+JuT)>$)~!@O z2fHAwJkY3_sW@JRai_@IwOopfgiSDVc;{eG1-DojxRchrkaF*c>m~2dsLi|QLGHDH zcVQWKqq1u@lzQ(!KuxfL+Fx;;QXI#q_GkZ1nd2ruyrfiDzVi`Yl@@#1*3Y2i+EYj? z)ni+`b9Jvm-G`~xr;(L@?;@Y+!Jk#okXDUt{eEH?B8H#o7~lkW&f0FT%G6?1^L~Yf z+1iJKlf2;8liY)4+xooy?RFb=xTg-0?+5>>tP?)pCcQ{`XP>VGQpWgzKxXwdaozJFP_)|>hR z`t3#Nfs$(uRB>89$SJB@Hg!M6z~LPP?xGH(z6=EsI}JnLeIPuG`lMsntvqs2eOal7 zHkpZTFbnY}>I~;lwtMG0%FZ{aG7jnvmGvn9X(Xx9sPq0zalSm>c|S*t)SI=w(eoMC zJKTm@e04_uL>q{{iW0~$biy(;94*;A7be2Ei^}fTNfqyy+!pc!cijmTOv9-1X*gdq zKpQ;=(J{*v=LuPOau+(~p(#7%i%^J8NqRTXY{vIl)}*~msm@aWYdxt{m#;?|TPsjn zZ3|-PjLPn>p~e1m2dTrYG<1-sLDF@<&FdiDe7I0~m=6@k5prbaMh(=aEFg!^-mja& z#`}?pX4zerc|U^i2~GZV<}CH2Stv*wPuE`BGo4o4a@YRU$WXe1lToa{sp*(sY5v_a zora6XLyNO+jpi$7YC3#t-uxk~G8+H1Kh96wi@=?lmZ_r!dE9!fJPfBQUglwKm|82R z5jVHdP`f>h>(QfUX-1)4kAX2w%hfy%wT^)I-!1;pWcn*B+cE^}UX@41|K*k0|fU1Z@?AG*r^eMzKMOcX-XBS*+ zt-T%3ipVZVz!d0Ec6~wYoE66vY?ntcd@rHkd+bfzEKRRh(_60cKfn7D z>}|1^$O+p)wP>-|q<0{j7lNWNKUJYghSLEKB1z?=Wzs7EU?X$xgS3_dcZ9HmznE3p z@3EIc;<@zoB*9HaxtLID@ZV&2;Qt<4_1aKe>pk`=&9qXuvF#`jG(NcrQxOWHZbjRS zc0z~J?wm%H@a{`4?wNmt0p2~e0T60L@>=<~$taJq@c^YFcB;PB90y40y163KZucsV z^Kc)~du#hhL6Z`F1;xUd@IcD-2MRkX1-rczjsnmYdl}7K$W0$fRgJa4K)bkEgXFU2HY?O9kHzuwSGH8Qzs(iOBmDI zZg=YpZ5=mM`IPw}?xpA1o1|%%?^o8VR4v-X4!WpIYtsgL;(FQ%U8H?bCZtX2nCEl# z6Z?ryGIo4Rf%;icKh$h0~=P_%92o2tG;MVU${OxGZb=5NoZn)>Tz z--PStPd27r@vLi1z2|B1q~2fW?Vn!;cx}_X%jDAZf1Sse8M5To-wK`Ykk-H_RJX+?N3-r3<0pfr1AL9w>OA;DLe% z3LYqUpx}Xm2MQi2c%a~czs3Vacy}?;5%(uMI=lS=X9w#DgnT{Ga3~h<^40iUtTDPZ z7*7Okv5;+U!sbuf=7wzH#4?Z>oHucIk3Sv^a9)y)LeL-K$%c457Pp0?wuHZ{I}(iB zLOs!rWH=UOE#a=9Ezuo}CP@+5-x$R!u3P;PTSqLK#LJ@L=oVYDGibvLu3Lgh8*+Mr zwrDJA3&nb(0a7ZJARr!HDvfT&dT+t+8F=CrRlN^)RPMVM>FSfp zVfUuxUoStBZM>V=la1u57jC{U^woX81wVxgf{xyD>06$!KXUx7X9N{E8{}+wWU+7m z*(;y${&4BC2Ww70w1ZxYV7(%b6Wm9+&yWT*jup)|-TW{Z=v~Tvzfe_r?$i5ARiy`5 z1EBO??WWgNbvy2cHlJ42$AP=wP*u9yvtLN2*=}9 zw)d2(Ruf14blgZ{%#Pn8{6=0=)lpDpv97dK-KxyoR@%#Mxq8_R3+>kvO!DpceHL}) zf&?NEp2V*g_y_?a5YL17RiW;9uFeBt`|z^^yOATB{0|X#16x3t*J2G?O1Vu(k?aKi z4b+i&EY`l_RhHgjQ(w^xi`8Q(t)E%^UCU>S8!ZFHeolMA_c5@Rkni49)sK-}gT?wl zvBy%?XSx-PrcDs8f-cCfoEf+5h4?;*6CEISvaz1z(H3}*fCs`_ft>@!eTOt*5n!cn ztLj3|g}Q?F0$T)ZKVe>rt@uA@SV|jaGTh>0UqN0PdA%h_z^=aHPgou=HZ__oPZl>> z24TjN#kT_4sMXIivp6MkIzc_i(dzdm^o*jt&XBb``#j(O-TW+aZZm}+>3Zrfiu|sb3tf~ePC;5Ly<*ObtHCSx@B@GtGgXVfm z%>&cCmL+}DS6SRq%M$2=Ug*!M(&`UArRP+YUV|!xf(Hs7D0raYfr1AL9{Ahwfc(x& ze%FPq1rJ6_-xbevgo_)V9uXD8H3m;hyYMKNr@tZc z@#Y>+%X*07{@+OMdu@Zp8rtw8xCLz%bd#Xng7ykJAm~0p9~X2`&=EmL1-&4sR?46M zFTXonwPM9Gn`7PPo@lbiw$!=ESu@|&!&9z%7t}aw7CBw>H2hENrY|w8CULL{7h&@I zbiEc6F23aV=lOUEyJu28=3M{I=bx7A*ZKH#Ci{OrUdm*@$j7f@vOneHGnnkR`S?sG z`*%K$(LPDLm@w*Xlj8UiagvV3r}{d@c-<97ike&a+)ZXvwJFhQ;qxLgmA5|>W@ri!>0&vnV*B7n?e_!SN6a7Dx{XgOSv)IG6Qo}f%tiI3l?8*u4_!`$U zu^rzM{Ko!pk>eBl)l5_zv%`dTr*N-XUBo(k-^k?`V)qa(7W&y zsx`t+Ou>If@Q(=o{eu7U6#O5BAIQ(6qCej#_-_JE^)-&irQC0VNk_a~j}8#g4Ztgk zKf>gCBRM)bzm*M&a^<})y6;y(Cn(ftoOH{F!tN>be_80UU8nH{HNyUV3jY5^z7^xK zs#!1!vOwsO=io>`|5V^^;ZM0QISici82iJ|1^+&=LF4zVP{uC=F6V824+{8cj_2D# zsrP}CACA*RgdqK(PQh=*ddbI&uv1Rov>wBo6T%I^sa=frF9vR#z>RJI{Q61o$Q1f% z{U&>?V!mAxj`~2i#|OE~#1Unw*Kws|fr2@wk77504|_J6I_0?+W?? zJzZTpAfik8z?sZr_4!ub;aT0_Yq+hR4&lgW2ryrL+ijlJjVtnV_?aDmaN>)v!7DU) z>+fK`RZVM_dzyS}R<2yz(Bf0sU8*aBeRD%x6z_}jdzyh47 z@yCu3p#we)+`)i9>DSKq$j{NW3mK4iCdfGFA{O^Y=}AyBh6hJDZe%=b!BRHKQ;G>k zgy;ldJv~m+KVCvB!#?dq5txy80!oCRaiWu72#!D^=Iew>I55TdSj*SA28sgVsIMmx z3~1+~=%mSa7l)qY^K3^~>kNmI3m1TWb0Q&Z3rG2LC5YhE6QAyAm2qimuI5HM zg=HM07873LXRzew_Jx&^!!FY0q4%*xk#&JpEfPg&98|NqeTNr2W^d~x*Gq59y zA~c$eYdQ4njvfheHj(flFCL8eNkGK9BT44u?chYzxg~}qMH0ad=1c~+BS|9}+%evX zoxx5qb~^*mCR3UYZ4he`GLEABU12EG#L<`xQ&3kICIvZF{WX-PVJx}W+h{>v7fC9w z-$gd(m2skfhaX)kOZg?3Pbo-hdr25U`dHep zp^Wt-LOzx9^7>QKLBS}?m-4dyj|lnYf=^zjN-FQ$kUYs4~VvitCUxAd9=raa+d0j54y#FWpjpOH_keBT*_rsF*<}x%5#`0e<$jj?|NgtL43u-L? zq(NTZ50KPYzU)6zEnOedHMf{{a`z47xqL=`nU?g|kR=yOd3k@Kxeg$LmqVY#T%1J+d4;b~|octQ%!t-qLM#OitWbR+eOG+mI=;d?QQEp8{ zmLnBNIhm&Kj0ltE%YC}lt#KInsVse{5arkCLByrJyv~!^QeICDU5vB3-I|137(mx; sbgd-qk@FAT&y}`o654;SSEG9jnX+8TOT*m2!n0;JX&hG@1Pn>`f14O`kN^Mx literal 17448 zcmeHP3vgW3c|Ouwwrs536|sVWfG;jC2I7@$jIjaNlI*pwK*o--z);|Ny$@@LeTdz= zMkYMO$V%7>1v7oblc61vS83{*xRarAW^h%A`NZ!D;)WY7q*V#JZ&n$Q zRCIX3l>Z{LZW0>E$vqYFcm>d zlkDm!t9UQ1n{VkcWw~9RX32$z?pMJub=pze>{w>m?O5xkhgG|j6lFVN$*x=3bt^lj zg7p>U_M|ZAZ?$SqFV!HUq{&|)ymmX3otJJ<`b?<|O1Hmf#{MhiZB=$+vE|20_bHeu zm-iLeQM>%MiMx55Dz92UJg3@`DfKs{=2UXms+G;D*s@eIQ|MjRyL#2KRVzdJY-oj) zQ(ROAjj2tyby2sJ2vf}{`(lCON&d6m#jE{&&zt@Tz+x{*JjPL)hRrHrbZp;j#NxU5 z?quGK=eBQNpUP(9+atSDabXy{)7gxXHzPUIFf7$+6;s*WBL6_nOo(_cA@a!#6%e%w zi1#K<(Ni$B%KCC}UG+NBeqXMQEj+TM@c`D!+y*NF&sTf3xIoof^SD;53%`dZ+=v3p|b67W^$p zCL3tR`l-SSmn7$s_6hKmNWTV@$9eJTKf&;?6TeWHrC)37^Zd&5TAvyOJJtNjc&oze zTj2J8+sqIF%f@x|K@N2?vf2x6?@n z?zhpI%uEZ+w7^UY%(TEv3(T~@|3(X(30(fE-haAYAFTVQW+C)}l380hsrSEJe@qr% zx&9Ww-`((fk8_~^e>DwX|`rv=OS=o>aq3#akXqc%?qm+{fRvUystjE_ER z^R!SIAAQ>9X@N36`nb*0!eo5(0h^}<$@pl*=4nATKKd1#r-j(~XdCm0zClL6LCE;Q z4}_>TQ~VqB; z6>6vtHN1$5>qBdY(O#9N3t@70KiDbleJzozKhuk2`q2;9=|?Agdd;i)8=sjA;NTVIpuRGaknP6y)c<{JNh}K2cIo|V zzfD=a_?|gmA6)x=Fy+0URVw8eI^)&4KLuZNFKq4d$$s=f2r1JAyEB2t@MI$J>r$XO zJh^8!*p1C#iLx)}Wdeg>tIKKRGv zG7-p>Dt`wLL*zKnd?N78QsA#k#erF9&cOjPcWDNiWdy8M53T+qq$|v({v*T+bG}s| zUC^g)R?@;3QHF|c2O_p2J$0F~DJ&o{m3Agzi7DTILZA$@2KXe3cu$6{e2hxd$4KVFz;ty-jruLU~H(9;MSbut4D? z4}>KAP$MZi0CTYepvpe9&;u0a9F}ceCo5GPSO{7 z$Udb&zXT)ex{=vE7jjW{j~Rq0wb$q>a;Z>!r*4)WLQ?CYW(Q-|p3=MyoxIiOvW*MD z$B$4G!{55VH=Y6MKOzO@HDqB8^TvTIp<;C)^@_X8PTg*ggNI2t)NqfS>Y9##xW~a) zKV>)6P@gvr^uf&&G5xyv^71q5%2&YBak&;BDeHsb@}MY&CkMlmM@Jf#TbXjB%8b+@ zBXznai$P1-HEGoA#YZN9Phgamzk~V}OPFBX1yH{S0bsnO@YM9@pt z1ruzz6Vo}3H%B!6r(59cO}$t+35O@i;UCHleza5!pC0Tw4fIXS056fG_#YVa#qgM3 z?4shw^uewPmoE*zWbH`%R@={6D3=LZZ9@|@<((~G|L@i`_X^s3JYY?2wXBSRfoI6WLPRMlXiE7M`$j@z>~9fc+wHf z!ZJz(mYzuHvaQH^0?P^c_`Te0+0;qhp|w9j@a45uBxv!jHxbFNpo~a*Wh5t6B*a;f z43RGsFUs%y8%MH6Me++3NyLhTj69KK>`1-`YibjzJG3@Wk^BfhXp<2{Qbr_c8A-Tr zFg$`V7NhI>tQh9nF&K6X&s+Y&@2VJnNyf171n>b?4pviKNv5g$EPIME#0kxvLG6rOL7}819@=$IE|%drBgJD{LSIS-;dkFr=u0 z`TO!)AIoYSg3TcW^PC;b58#=DISlDx4rZ1e%!{yMg{TnLwc+IDB>(Au9} zAw=699mRiW>*^?e+_t@~IN7BSEqej*?VVSjqcwVY%UP6t^qkLJwEQi#e(NlbcNRa~ zQ2bR}Wx?Be|FIf<%`Xe@)B64Q?`^xM?cTQg+Kgj~M95hlt?a2>x$d`D`ZAoEpP3ez zX@UPoEkN&CDCMKMh#BqfiNr!t;yW^UMYSiA(xTaniFZQD%x=x>j%#>9wL5NVkQCxt zCTnVmY#|fVk{K;g$VAO#HiP`0crG6|qOru6Nw6fZEnO13n!I8zi$`+NZmlPmjmGnN zEtZT)FDL~@yW`T}xAF8=praK97==tSkn%jK}7o zS~4~*n;>b#)RrW)Wd25wd4k;KA1Fj}@tAZ+d5{LvzSS7~EexpD5RBeZU;WvC?zcw^0>|gq1-3EWl*JkVf)%|le z`CBvo)ouQkHviIgzt-+w*zOOu`|IU<-Lvr113%99l;m3l-m_C=+ESy48+EE$cewWVFV3K_GYtqH9RwJcj+klgY|SG0s$R)&^eZQ;MsCBEjYvA^$kYSXu&{C&C` z_X+;~+>Or?4^PuEK8#`>&u;zMg2%BNpCfp@yYYI#>0CLhh#B zHa^iLcs_DUeBz6O=OH(Kk(ll}@QH=?yy4cztk>#9th`TLB7)Q6JS$Gw(DXA)d{mi2 ztQ^g^wdjYNM9@$OyXLVYL9_!SE0`JC|^kg0LRL*-H0h`*X?+93^`%I)*`-z)Xc zZ_fk3r;EQ@y}l>?ED%r90}@Y~&d-mePSbh)^)m1paejZjuJpbAJ|gk+$Ilt*=luT7 zVQPsb2Vz8c69^YEK3*Sg(~_#SNd-${UN3Q zb`|~4l>TdKe9C8Bu<@fkseQcH+YSF*iT^k0NLi0gD=tHTMh3zQ!|@7AouLkh1^?bfF71Hj3T zcbt4z>36FIihOnko4)`~^*XFs4!B+irM`0NAV-v+r&PI}O7=<>{eL9?XlHLc?*P{v zemQ#%_+n-p^kUp!sb3n@cn-5Hnltle0dEkaV(RRQVWtfnDapj?s7Wkq>`rBO;W2s4 z%;xe&q|hs(S$ZyxNA983EvuGKuS6$Ml13z#i|jSAf&$bMd9qZi^dDaFHPF@*dU3t(qY!k+&&Mob2 zoyL}p8*dM9H@3I6cZSKUJS7v&7i8RY{HC>)&&SZIn9n=jLx*xa+VNP#j9BM#+!Dv@ z6BOz`$1|n5(mtd!g;5l_NQT~Cm|46J!C4`PaE+MeeZqNXe;fwz=6xlC&PLSMQ6U89 zKC@)R@>!!B<>Ne)_tlKiu?1GKWX34u<1y<5le;KK`%-7Cvr2BAUJS49WWu~|1q$Dl z&#OS>nWe2AbQ;R7uPQ9)V3b>87&v~#Z{{r8 z%{d&UisyPn2rr}Jq1~B62yef7;yH7#1KL$c;;dmZrl7X=j%8+Kw~*@Hk$krZ#r9@U zghkDqC83vDc;)6|49IfvRD=vvt|w)RkZjcuvZ38su$0NiqatL+d%@DR2<@yKm7#dI zy860f@W!0w!@3MD3(ljcNID5;mN~jcgfIo9aar>`^53L1&EwXDHG`EFs{H^{-tQ{0 z)a5+U-@za4G1;E?A51mKoQe09jM2G2r~M99KGQ2%QBj`f8Se1d^L~b@#)^u1?e7GR z*JoLdpZ7gXxqP*qwfFt(Hi8W8;n|+|M@$b;<04TVDa-Ze{dN)|P?8(M!|128ebV~PY0vvwrVXs9sJH!I zQTE*bLQ{%NZGnYz|0yl;Ix~i_MkLUlVfE$LcKT2%J{rH!_ zsQ$M7I?}>70cyd0z+&xRwr4sU#!h?Ds#IGn%H^;Fwqt%ij0xlNdA%G|_FjD|%b6Ng z`AXduvHi6+VrAIgNj+UG6 // +#define STOP 1.0e-8 +#define TINY 1.0e-30 + +struct box incbeta(float a, float b, float x) { + // Descended from , + // but modified to return a box struct and floats instead of doubles. + // [x] to do: add attribution in README + // Original code under this license: + /* + * zlib License + * + * Regularized Incomplete Beta Function + * + * Copyright (c) 2016, 2017 Lewis Van Winkle + * http://CodePlea.com + * + * This software is provided 'as-is', without any express or implied + * warranty. In no event will the authors be held liable for any damages + * arising from the use of this software. + * + * Permission is granted to anyone to use this software for any purpose, + * including commercial applications, and to alter it and redistribute it + * freely, subject to the following restrictions: + * + * 1. The origin of this software must not be misrepresented; you must not + * claim that you wrote the original software. If you use this software + * in a product, an acknowledgement in the product documentation would be + * appreciated but is not required. + * 2. Altered source versions must be plainly marked as such, and must not be + * misrepresented as being the original software. + * 3. This notice may not be removed or altered from any source distribution. + */ + struct box result; + + if (x < 0.0 || x > 1.0){ + if(EXIT_ON_ERROR){ + printf("x out of bounds, 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__); + result.empty = 1; + result.error_msg = error_msg; + return result; + } + } + + /*The continued fraction converges nicely for x < (a+1)/(a+b+2)*/ + if (x > (a+1.0)/(a+b+2.0)) { + struct box symmetric_incbeta = incbeta(b,a,1.0-x); + if(symmetric_incbeta.empty){ + return symmetric_incbeta; // propagate error + }else{ + result.empty = 0; + result.content = 1-symmetric_incbeta.content; + return result; + } + } + + /*Find the first part before the continued fraction.*/ + const float lbeta_ab = lgamma(a)+lgamma(b)-lgamma(a+b); + const float front = exp(log(x)*a+log(1.0-x)*b-lbeta_ab) / a; + + /*Use Lentz's algorithm to evaluate the continued fraction.*/ + float f = 1.0, c = 1.0, d = 0.0; + + int i, m; + for (i = 0; i <= 200; ++i) { + m = i/2; + + float numerator; + if (i == 0) { + numerator = 1.0; /*First numerator is 1.0.*/ + } else if (i % 2 == 0) { + numerator = (m*(b-m)*x)/((a+2.0*m-1.0)*(a+2.0*m)); /*Even term.*/ + } else { + numerator = -((a+m)*(a+b+m)*x)/((a+2.0*m)*(a+2.0*m+1)); /*Odd term.*/ + } + + /*Do an iteration of Lentz's algorithm.*/ + d = 1.0 + numerator * d; + if (fabs(d) < TINY) d = TINY; + d = 1.0 / d; + + c = 1.0 + numerator / c; + if (fabs(c) < TINY) c = TINY; + + const float cd = c*d; + f *= cd; + + /*Check for stop.*/ + if (fabs(1.0-cd) < STOP) { + result.content = front * (f-1.0); + result.empty = 0; + return result; + } + } + + if(EXIT_ON_ERROR){ + printf("More loops needed, did not converge, in function incbeta, in %s (%d)", __FILE__, __LINE__); + exit(1); + }else{ + char error_msg[200]; + snprintf(error_msg, 200, "More loops needed, did not converge, in function incbeta, in %s (%d)", __FILE__, __LINE__); + result.empty = 1; + result.error_msg = error_msg; + return result; + } +} + +struct box cdf_beta(float x){ + float successes = 1, failures = (2023-1945); + return incbeta(successes, failures, x); +} + +float cdf_dangerous_beta(float x){ + // So the thing is, this works + // 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; + } +} +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() { @@ -215,7 +347,7 @@ int main() // set randomness seed uint32_t* seed = malloc(sizeof(uint32_t)); *seed = 1000; // xorshift can't start with 0 - int n = 1000000; + int n = 100; printf("\n\nGetting some samples from %s:\n", name_3); clock_t begin = clock(); @@ -224,11 +356,11 @@ int main() if (sample.empty) { printf("Error in sampler function"); } else { - // printf("%f\n", sample.content); + printf("%f\n", sample.content); } } clock_t end = clock(); - double time_spent = (double)(end - begin) / CLOCKS_PER_SEC; + float time_spent = (float)(end - begin) / CLOCKS_PER_SEC; printf("Time spent: %f", time_spent); // Get some normal samples using the previous method. @@ -236,10 +368,10 @@ int main() printf("\n\nGetting some samples from sampler_normal_0_1\n"); for (int i = 0; i < n; i++) { float normal_sample = sampler_normal_0_1(seed); - // printf("%f\n", normal_sample); + printf("%f\n", normal_sample); } clock_t end_2 = clock(); - double time_spent_2 = (double)(end_2 - begin_2) / CLOCKS_PER_SEC; + float time_spent_2 = (float)(end_2 - begin_2) / CLOCKS_PER_SEC; printf("Time spent: %f", time_spent_2); return 0;