From 306748dcee4955d4dbcebc5a61c47803d64bb61e Mon Sep 17 00:00:00 2001 From: John Maddock Date: Sat, 2 Dec 2006 19:21:13 +0000 Subject: [PATCH] Small fixes to beta and beta distribution. Brought docs into line. [SVN r3487] --- doc/distributions/beta.qbk | 129 +- doc/equations/beta_dist_kurtosis.png | Bin 0 -> 853 bytes doc/equations/dist_reference.xml | 71 + doc/graphs/beta_dist.png | Bin 0 -> 6265 bytes doc/graphs/beta_dist.ps | 12004 ++++++++++++++++ doc/graphs/beta_dist.rgd | 1539 ++ include/boost/math/distributions/beta.hpp | 41 +- include/boost/math/special_functions/beta.hpp | 13 +- test/test_beta_dist.cpp | 219 +- 9 files changed, 13831 insertions(+), 185 deletions(-) create mode 100644 doc/equations/beta_dist_kurtosis.png create mode 100644 doc/graphs/beta_dist.png create mode 100644 doc/graphs/beta_dist.ps create mode 100644 doc/graphs/beta_dist.rgd diff --git a/doc/distributions/beta.qbk b/doc/distributions/beta.qbk index 3e1af3a6c..29dcbb307 100644 --- a/doc/distributions/beta.qbk +++ b/doc/distributions/beta.qbk @@ -45,18 +45,6 @@ RealType probability); // probability cdf. }; - RealType mean(const beta_distribution& dist); - // Mean of beta distribution = a/(a+b). - RealType variance(const beta_distribution& dist); - // Variance of beta distribution = a * b / (a+b)^2 * (a + b + 1). - RealType mode(const beta_distribution& dist); - // Mode of beta distribution = (a-1) / (a + b + 2) - RealType skewness(const beta_distribution& dist); - // Skewness of beta distribution = 2 (b-a) sqrt(a+b+1)/(a+b+2) * sqrt(a * b) - RealType kurtosis(const beta_distribution& dist); - RealType kurtosis_excess(const beta_distribution& dist); - // See below. - }} // namespaces The class type `beta_distribution` represents a @@ -72,16 +60,6 @@ See also: [@http://documents.wolfram.com/calculationcenter/v2/Functions/ListsMatrices/Statistics/BetaDistribution.html beta distribution] and [@http://en.wikipedia.org/wiki/Bayesian_statistics Bayesian statistics]. -[h4 Member Functions] - - beta_distribution(RealType alpha, RealType beta); - -Constructor example: beta_distribution<> mybeta22(2., 5); - -constructs a double beta distribution with two unequal shape parameters, -whose pdf is shown in grey in the -[@http://en.wikipedia.org/wiki/Beta_distribution plot]. - The [@http://en.wikipedia.org/wiki/Probability_density_function probability density function PDF] for the [@http://en.wikipedia.org/wiki/Beta_distribution beta distribution] defined on the interval \[0,1\] is given by: @@ -90,14 +68,15 @@ f(x;[alpha],[beta]) = x[super[alpha] - 1] (1 - x)[super[beta] -1] / B([alpha], [ where B([alpha], [beta]) is the [@http://en.wikipedia.org/wiki/Beta_function beta function], -implemented in the Math Toolkit as __beta. Division by the beta function +implemented in this library as __beta. Division by the beta function ensures that the pdf is normalized to the range zero to unity. -[@http://en.wikipedia.org/wiki/Image:Beta_distribution_pdf.png Examples of pdf for various values of shape parameters are shown here.] -x is a [@http://en.wikipedia.org/wiki/Random_variate random variate] -in the interval zero to unity. -Note the [alpha] = [beta] = 2 (purple line) is dome-shaped, -and might be approximated by a symmetrical trianglar distribution. +The following graph illustrates examples of the pdf for various values +of the shape parameters. Note the [alpha] = [beta] = 2 (blue line) +is dome-shaped, and might be approximated by a symmetrical trianglar +distribution. + +[$../graphs/beta_dist.png] If [alpha] = [beta] = 1, then it is a __space [@http://en.wikipedia.org/wiki/Uniform_distribution_%28continuous%29 uniform distribution], @@ -107,8 +86,22 @@ If [alpha] != [beta], then the shape is asymmetric and could be approximated by a triangle whose apex is away from the centre (where x = half). -Accessor member functions simply retrieve the shape parameters, -/alpha/ and /beta/ used to construct the distribution. +[h4 Member Functions] + +[h5 Constructor] + + beta_distribution(RealType alpha, RealType beta); + +Constructs a beta distribution with shape parameters /alpha/ and /beta/. + +For example: + + beta_distribution<> mybeta(2, 5); + +Constructs a the beta distribution with alpha=2 and beta=5 (shown in yellow +in the graph above). + +[h5 Parameter Accessors] RealType alpha() const; @@ -118,34 +111,65 @@ Returns the parameter /alpha/ from which this distribution was constructed. Returns the parameter /beta/ from which this distribution was constructed. -So +So for example: - mybeta25.alpha() == 2.; - mybeta25.beta() == 5.; + beta_distribution<> mybeta(2, 5); + assert(mybeta.alpha() == 2.); // mybeta.alpha() returns 2 + assert(mybeta.beta() == 5.); // mybeta.beta() returns 5 [h4 Parameter Estimators] Two pairs of parameter estimators are provided. -One estimates either [alpha] __space or [beta] __space from presumed-known mean and variance. +One estimates either [alpha] __space or [beta] __space +from presumed-known mean and variance. -The other pair estimates either [alpha] __space or [beta] __space from the cdf and x. +The other pair estimates either [alpha] __space or [beta] __space from +the cdf and x. -(It is also possible to estimate [alpha] __space and [beta] __space from 'known' mode & quantile. -For example, calculators are provided -[@http://www.ausvet.com.au/pprev/content.php?page=PPscript Pooled Prevalance Calculator] and +It is also possible to estimate [alpha] __space and [beta] __space from +'known' mode & quantile. For example, calculators are provided by the +[@http://www.ausvet.com.au/pprev/content.php?page=PPscript +Pooled Prevalance Calculator] and [@http://www.epi.ucdavis.edu/diagnostictests/betabuster.html Beta Buster] -but this is NOT (yet) implemented here). +but this is not yet implemented here. + static RealType estimate_alpha( + RealType mean, // Expected value of mean. + RealType variance); // Expected value of variance. + +Returns the unique value of [alpha][space] that corresponds to a +beta distribution with mean /mean/ and variance /variance/. + + static RealType estimate_beta( + RealType mean, // Expected value of mean. + RealType variance); // Expected value of variance. + +Returns the unique value of [beta][space] that corresponds to a +beta distribution with mean /mean/ and variance /variance/. + + static RealType estimate_alpha( + RealType beta, // from beta. + RealType x, // x. + RealType probability); // probability cdf + +Returns the value of [alpha][space] that gives: +`cdf(beta_distribution(alpha, beta), x) == probability`. + + static RealType estimate_beta( + RealType alpha, // alpha. + RealType x, // probability x. + RealType probability); // probability cdf. +Returns the value of [beta][space] that gives: +`cdf(beta_distribution(alpha, beta), x) == probability`. + [h4 Non-member Accessor Functions] All the [link math_toolkit.dist.dist_ref.nmp usual non-member accessor functions] that are generic to all distributions are supported: __usual_accessors. -The formulae for calculating these are shown as comments above -(and in the table below) with the exception of the __kurtosis and __kurtosis_excess -whose formula is difficult to show as a comment but is given below, and at +The formulae for calculating these are shown in the table below, and at [@http://mathworld.wolfram.com/BetaDistribution.html Wolfram Mathworld]. [h4 Applications] @@ -185,22 +209,23 @@ please refer to these functions for information on accuracy. [h4 Implementation] -In the following table /p/ is the probablity and /q = 1-p/. +In the following table /a/ and /b/ are the parameters [alpha][space] and [beta], +/x/ is the random variable, /p/ is the probablity and /q = 1-p/. [table [[Function][Implementation Notes]] -[[pdf][`pow(x, (a-1)) * pow((1 - x), (b-1))/ beta(a, b)`]] -[[cdf][inverse incomplete beta function `ibeta(a, b, x)`]] -[[cdf complement][`1 - ibeta(a, b, x)`]] -[[quantile][inverse incomplete beta function `ibeta_inv(a, b, p)`]] -[[quantile from the complement][inverse incomplete beta function complement `ibetac_inv(a, b, q)`]] +[[pdf][f(x;[alpha],[beta]) = x[super[alpha] - 1] (1 - x)[super[beta] -1] / B([alpha], [beta])\n\n + Implemented using __ibeta_derivative(a, b, x).]] +[[cdf][Using the incomplete beta function __ibeta(a, b, x)]] +[[cdf complement][__ibetac(a, b, x)]] +[[quantile][Using the inverse incomplete beta function __ibeta_inv(a, b, p)]] +[[quantile from the complement][__ibetac_inv(a, b, q)]] [[mean][`a/(a+b)`]] [[variance][`a * b / (a+b)^2 * (a + b + 1)`]] [[mode][`(a-1) / (a + b + 2)`]] [[skewness][`2 (b-a) sqrt(a+b+1)/(a+b+2) * sqrt(a * b)`]] -[[kurtosis][`3 + 6 * a * a * a - a * a * (2 * b -1) + b * b * (b + 1) - 2 * a * b * (b + 2) / - a * b * (a + b + 2) * (a + b + 3)`]] -[[kurtosis excess][`kurtosis - 3`]] +[[kurtosis excess][ [$../equations/beta_dist_kurtosis.png] ]] +[[kurtosis][`kurtosis + 3`]] [[parameter estimation][ ]] [[alpha\n from mean and variance][`mean * (( (mean * (1 - mean)) / variance)- 1)`]] @@ -209,7 +234,7 @@ from mean and variance][`(1 - mean) * (((mean * (1 - mean)) /variance)-1)`]] [[The member functions `estimate_alpha` and `estimate_beta`\n from cdf and probability x\n and *either* `alpha` or `beta`][Implemented in terms of the inverse incomplete beta functions\n -__ibetac_inv, __ibeta_inv, and __ibetac_invb respectively.]] +__ibeta_inva, and __ibeta_invb respectively.]] [[`estimate_alpha`][`ibeta_inva(beta, x, probability)`]] [[`estimate_beta`][`ibeta_invb(alpha, x, probability)`]] ] diff --git a/doc/equations/beta_dist_kurtosis.png b/doc/equations/beta_dist_kurtosis.png new file mode 100644 index 0000000000000000000000000000000000000000..1be258bfdeca74a4e9173205c21d1d34e1d97e12 GIT binary patch literal 853 zcmeAS@N?(olHy`uVBq!ia0y~yVDw~QV9@4ZW?*1AB;Qrfz`($g?&#~tz_9*=IcwKT z1_lPn64!{5;QX|b^2DN42FH~Aq*MjZ+{E{6KPd0WZ5Gp3rwv7{Pu{Y1>wA?5Lz$mXe#Wvdy8U8O@TKdl ztMbpj_6fh;;4__ZLviGTExQ-2xGi%;_SVJy>*Z`#JnmW5c5=sW#xp9L?)~vqy|j9j zrLJ4tmIr@k-gRcR>R-=%-B&L3E$^{U=9gY?|8j2i>}~{rsd9EgR_7m1UVPqGPh>a;Qt!f64fh zo1`Tt#3tN7qw;CG;mwdYr#{IK|y(Te_IqSglNGEvh*yt + + + 6 + + + + α + 3 + + − + + α + 2 + + + + 2 + β + − + 1 + + + + + + β + 2 + + + + β + − + 1 + + + − + 2 + α + β + + + β + + + 2 + + + + + α + β + + + α + + + β + + + 2 + + + + + α + + + β + + + 3 + + + + + + \ No newline at end of file diff --git a/doc/graphs/beta_dist.png b/doc/graphs/beta_dist.png new file mode 100644 index 0000000000000000000000000000000000000000..3ddb14f9ff9412270cff9ae956a2bde8c8c4e9e6 GIT binary patch literal 6265 zcmeAS@N?(olHy`uVBq!ia0y~yU@~W5V4TRo%)r2~E;{KG0|Ns~x}&cn1H<|g=B!;W z85kHOOI#yLg7ec#$`gxH85~pclTsBta}(23gHjVyDhp4h+AuIMNCfzVxNhFG>B0T` zc6N4auCDt3f8(lE_y2EL|NsC0#`hB@Gcbrh_jGX#skrs#ZfDctRYGhJv`l{e%fG{T zCuujsCXpZ4=Uy&fW}O?dPo$&4L`mSEqXjF+4#xw$OnCwY%??nyL7L@`f~duPkzwtXMDpVtST zU;pc0gPD5kH1WFq)eqACSMA6=(aY*!o_FlO{_W@Un9cRqzKDN(?4SJ2RUS9@7=`;p zD=f46b%A%Z@qU}Z`;X)tW7sc4$oedUgrGjIcv)2 zOLID|Y-?O|<&*sR`w?+W=4-cfmnJE=dUstAch1iGF+Kkc{~GyQg*SIBd>LNxUwB96 zvRY^FSAoUC&rY1Zv@h`Sl>E9K7aaE~n!c-NGC#W`a{bZeYA63!buF(k-nwvN{sT8j zk?h}Tzu4c#rQR}`{q0o$?A+gN{kP2vF8lp<)4hIw+0x|V-hCmrt_9By@4d|bwtHLN zGHrw3;;D20#`z!gjlLed=g-|3H{ZGXdv+Q2&i&bx{kH$2!ijI!KImlcQ!?Mb?frhq z`xm4;e$N6&=b}RMJ<^LC8l=_D*%wC&6g-|#;8A{4Ijc|TXWBikg?m*L*}l049O_lt zA^0?H&7ret$+CB*@IU)?XF+RB;x5hiX>T5|zn7l*tY1*30zNZTwsBN8X_Aul@ zn#Aj}%PF@d@-FSpG<Otp41r&(H)$(^fs%Qzde|FmA7}%9h2+5%J{GNJRJJMG_=-qqy8(PzzFWcCyRh>~&d;Q`WMfs~+KHdK=zh-~3 znEC>JHtS8j+g@||akFolxh8M=mJPx@&saI<-tAetBa^%D_RdS)6`AIL{w}`0?(eg{ z(v4FzuhvY=-M2LU%qrvcOMH2kAItn2z5V<1wmT-Q$6v<0WYy`Bj(t{jk^Aq|DY5fQ z$_{6Trrw>Aa^&*+S*3<*=6hDI((v6S_;T6O)rIcni%XribtbKvW4ZIi{#U|_uH|Ry z78fOLRhgcbb$>#XwDbC+%K^sr!2E>R+{_l*J|5(WoG9u zGdEuS`p=rJzL_^TN-z5z`)C#ZT|2Mv&Zew$zQzSfyC2SYfBw##*(}HJ}C?7uiQtb8n|HanfX^vla-p{XMYuVYl;EV2xvyZK=mtLRw zJatFu`+^4Pyc@oM3zjKM_a*O?U-$pUZ=N;Peb$#VZITzJ-r@EAE`G+B{jc@zn5V&E z{Wmrgm|s0yc+UKEt#9@GgEl7qzP_{X{`INm-aOl=_h(|Bxd~T!aJ6~g>ssHqV}UFD z9vhZxm&zW^dhfQK{p(r3x#bm6*Dl{bYq0I@vLcJVyk~{c-@>-l7Vzj!559lwOpX0z z8(X)Tzh^F)cP04j?v=Y&70GV@x_61mW2G}kRcEi?lIi`yM0D2mZLR<2K0UMH{crx5 znbPy--QOD7{w=RN_H*95O>RG4eVU?jcS?aV@4hg(>VQH}NY8oucc+#4wyDe8kNPd% z8oTYr-ow88zRMqbb?1Xk#IB??Q^PYb> zDP#MbHt*S;`omonudXXrcP@OOX7YFrv!+J+yfc?d-Rt|F{k)VmZ~ecdJT0Bdqu<2W zf3AM}SJSSx?Z(CL>~<OIu90=sNjv?(^*BcRr@gTkmI_ZM*c0FUQ*!^Tp~OKcCE2 zOkc=ZP|W&V%uUhK@pDb4_4)G4GiKh=xb#?O+RweqzI`eo}*Yprs$X+w97H zd*(miU7CFC{)PP?%}=eaHSSqoT+DPl)P>P<-{+dMMP~M&Yuta&eY)%9x!-@c&AYx@ z8wWqOO$hJX@$XS$ z&C1#F7K~h*4mC;}=4~rD#FH{b`RrMPeb>)g^=XSZ*jW@PpE%2-zURm0gOzsYr@u_N zxb=RS(R_RU!tei>w?9?mIJfxF`=dK8cgp99b>4~3|GQ^J@N;1^P1jv3ulF{7Ui0Gf zk~!?3Yuu~nH=2E3{;Y2Q`}(_9=db^@WBz>mZp}j*-S0XD%SzY({A~H_vC;YRnEV{Q z*B$?6KYv{@KT#tpP3cb0-N#Ai8J{1tiT6F>!yCT;yXo@EeOqQeKb`maokZTg)_*l+ zK5Q2EKG&?3HqO7yck0%i2Y>9(U)Q_le=7IhU4wo1di6eWTwYbBn78eG`iyJd^V_e? zcvEC~e!9%2m$FA>wS^uk$ z%m*j5y{HJo??NY`1d5?FW z`CKY_-mWLVyg1}g(ZdUORIGa{zptP7-1^@g{`>aFuiu$p6cPU5YhK@b&CXVkzGpjU zdfRTbIbMEc&9nRa_R8D!tS>Jvde(ZV@}JH9uQSitpZ|FG%xBs6$M0`>eQe|7_ge~D z4|RU9xqnvc=IlN5NT+qnxPh*VIv!DC!o%P+g$;LwHS9kkHN%KAH4hmlMZDF+B z`0#*ljQRikRY5FfpH^ux=Cy5i zh;er9-Sh7YK5uxvqkVS6qpX8x4__M#65p@Gc|N>-lfeaAGdYys+6ANlEby*;_IC)Z$0NERg2FlyEEhFd-jsL z?=6?7-kRJuM=qyp=f&6RGByoO@62S^o{`(qr=2ZbtYyBxyvgWlaLtT>Pdh7F*k2x7 zx#7|7K=rEuZsz_{mZk?iOL^lC^hOrFDfoJ5cWh^L!uk^yColP~ta-e_&?_cR|AfWM z8wF+NNAiBE?VSDL^s1;^$L46?<#?t4Z~xmLtw(zgdmdlbAicS7vy#|@XN+^o&AYmH zEbLpx`eeQNws;0Z&&XD@qO69=26qeU%zxe7a_Z>WclSHJR^9EHJmIWw;L?-kWsdxs zc}~WCKh*flpQXjjD=yl;@&39uDNjWzB2{Jmyv$PX+)(GU$m#iXg}ZTW^3A*;v#$UC z!KcdSH$K^sDdcypWb*SWua$lw`#AHSe9=!0KGK^hy8FxhKQV8Y>Mc^RD4K2InZ2_4 z=9Ztox?}J3{Jhh$weRSmlXLQlwst@N9lUGlVfNP-#0rG3op^cw#Ih;Usq#BbZQHx2 zEib%S5fzxKIe)?~FBKWT*KIo|NdKO)?dA2a%qr!t-~7K{H}6p==kvq!`+xe^-~7Ar ze!j-tlv(VH=iBeT7W{wy>y0<<(*H)Vn-rD#&XH_w?()}i$do?%^iEFw+?zIU-j^g# z*>3hA!6M9iLZknW^3O{M)mhSkGx{NuRc6hWBlH=JNtEo?!z0uziWI{@i&h-#%p;lUE|rY zrKNvniw2eaFI+J9O3pcN&E(xT5|ejvt5>U^Og&|}!c*w0^yd1xQ};}JRFwB+uJ3oR zO)ayFE!0(K=JDHJ_M2HUd1dYAy??^q?LMOVE2aO{s@47VmlEaPl}y_Fa8jPz@sH-0 zw|9Se`hNZDcU<~!{j4r|#(B@37RYCQ?(tu*T|ZZ@j_K0<8fW!}S>4~<_3jI4UODmN z+Li3H{0e03UrQbfx%1YO%`E0Z@2A@l{xMvB<+lX%n2%r3jPRa2t$bP8#@jYCmNIF4 zV|nWLGlo5VmH&B_;+-HTJ+^qCr`hYo>3*Fn&RP4{oa(Rz)2#I0$>iNwKk3bJrT(1x z+AQlAiq^Vap3juCV}W#7yLj;Rydw2EXW4H(cT2za@!0+%Ba=KmFYD7aZHLdz_BO8+ zKDX}eNx9O{dyTJdzgm^C>*PbXpF3x|xdv*Q&MOUlwr^Q#y77jcTFpK`-G1&mIjx|t z+9N_|j;G(vD6dBrUH{~&vrN^VZV0icTHn5D>y#sfB70BgEzihcyB_-UVZo`Bf{RbH zhtHX9cQ*Rxs*fvIfAW7Mv(6~BEy`=l>@B8K5367Mxl7Aj=JpQH*~NL6_9mEXW@{ZU zeis$~v~=^Gc#j#y@9ef*y7NWuS@7zQD^~v!$_aWVBc)w7d0FS<@^!MWSM?OMvAYYe zoHps5nxpV>>+ROpLT^6!^LnSJ!t;N7PL~PqUG`{#n4@}%$aTZnEJZSQ>ki%da*xwE zQS|BbBBO(ni)v-$L^qX~8$N8Q*NNnrqxfrMujBHVJ+m|3|9Ep|*OgGoJ27{rb*UYH zIa%wz@0@KPAN#xu3jfM_(`EAGX%A9&zGRNh?3{cyB)l-wqMW%a*`!$JHmm54NzAIZTBb#Tp7g*do6t=w}Xj|{yW%@RgAJ>Qmwa#`gSX7pGXW5*Vg0?U6 z!mo}=-JLVPm~XF4|2@?(p1$84`yV8Kvp6|*2j|=4o&|@0=si_wdTeT7@?yhYyV+NF za5heUHt&AyDbGj8iu2MAO#Ho-)oVwh=%Z&jK`xh8Us^EbvUs}igd&^sd7S3yuit#% zKA*$y_r&h0MZZ4f)b!sspDn&}d4^3)`u_fkk~tMW-j<%~e8;sKP~($8`-BN8J2WSpEA2( zvA_GAFSDKeZcIFS-b_um{y|#m;fZ=eLZM-d?PAuqOFr-HEd8o*^lbC9v&%9MKizyT z@YF~99V>F`>;np;ny&7?nI5<7Qs#4qzULQ@zS;h>-fGd6(5<^I||J;WnO;V|LLnW^+msWa^Cg-E>w2! zn5|&uCu%y?%jWl~Qy&H7*DZKID}BCMkyhE`k25u!(&rb=>HFxoEph#qUmMn~y53v9 zYU!J+*ZlepJF=zCp5R=~wmizb(bVMm)z8|`ziI4cZ@qPxQ)i0**T>uT+_X8f>GI2S zcG1iFkE=g0lh0WQoizjDFeB5?HNh9j1%;A|g z{_M$|H?P@$Z{Dvj$0b7p12*av>CCw+#;$eEQSyuBnzOv(a&9#@o)z@u`NiFx=h1h3 z&9S|wR-EZoF;ko#b~iE9R75UctRU=l(dU)3&OTeTM|^S))5DOb%l<5Bm6#s(XLb?G z2F?cMgJ);o&?!o@XzpqL)EYc-w%d&}MLaTcM`Whlf4TM;@AkuKlRmd9FY4oZ^ zYjPqjOQ8PgiKfpJU$lAtw|jf-=KJ^)X}LEG?stCrBlzOk?DQ#dI{#m%Z$4M|-^uLu zuLbXS{dypB!hBCoYs`(BS;?FC_q=jFS=c;zZ|&hTUX?pCSKWSAziH-!4GNZi1*hWu zzrIf@zLWlW@48o3>p1!f%-3$T+`E2({)wX}3LBJf`)22?v+6keR!i~fRogj_^5Tv> zd$z6VE%Ow`=VE;+=I@g>^Q;vJjJ>P(Rq?uSbk3Z2iF!p=FAJXR%IxZWJ@@_7y3k6? z>jg)4Wp?(Sf3c=l=H%k5OF^9b2D874RlZ!kTH9jjU6b=|f7#P>*}|vzT8S5=?8d2XJv8yfGEuYsHva!79wvRgDk;;LcPrtr zKl|MGZ4;)gTJQHjr0lYmbe~~pqWzB@KDDJh?sH~~-6&YYsrSEa?erDv`5t-{UuKxC z_-xCH`XVlq(!K{A#liNoUlhOm=RftN-*bjO9rLZ5AF3D|@T*0pE-U`<ds^j;X zqJn%^ZL*m)Z~D}8`zK9&vMkHY`PRvM?ZVc^mu^n)RWKLb{7~_4m5Tc73D0&}+@IPf zV}5nhL&eLxgwiLU{CweOZgp2-u=Z@l=SeDmS!K_0x%bSLy7^&3J?retqB@tZM?_n= z-YqdrpDmLswSUQEo3}X+I4b{zu5f>yf_r zB+sVujr*hHjMA2-J`^ce{>in6$!J|>x|fX%w_m2Y(@h@Hi61v4e=HGvFr`enmp5GV z+4Xa8nbzKQw!Od1HPXCkQ}0(J6T|6?{MSA|oO-urNx@2s!!8#e-~RAZq%`g-hS2Gpsaz2f??L65UAef{otEP>2pxGQU$1mWE_m-e%g?F&t7jixdN*WQR^^v%qWg9nd)zF&qY;8;w@<&+yMM>A z$Ey$WF8b6U?V_8vY<5PT=akC7y*u>P`eX98dmcGgZCG1c{qax7$sfm$-|3S+dTNgL z^1sH*LrPD5Ji4cH&xf;dH*&4^MZ7oKxpU)gpFFWQW}SEY^S@Zy6^TYRO20b8pdZxt zd-#^t;lv1;3eHdi~I%AS7n@6jiPHB(l#c<)gBy)JLQ)NiX#RT~$3 zOg`yWbT75+cDkvk!NlWw{qK%1x9RJbjgM)TR+}$!$js7GuBT4Pba&U4sk5T@ZQO11 z+-`Yn*`Fn|W4~G3eX`g$?Njf{WwQ3WcUisO*ZOSQm;LkgYQyfv-O()izv6amnBUdT zyQ$a3XUC>(Ppo1*{(PJG4nFq<|3L%YI~>6S-JaI`Zifz}?cjW26JWtz@K}~B=K)*Z z1GQ+|)?!xk4DrK|0rLiF7VyA1bfCNuN`nW;!RCNRR+;hyU~~a!Or26KJNoPy*%*?~ XE<7d_w inline bool check_alpha(const char* function, const RealType& alpha, RealType* result) { - if(!(boost::math::isfinite)(alpha) || (alpha < 0)) + if(!(boost::math::isfinite)(alpha) || (alpha <= 0)) { *result = tools::domain_error( function, - "Alpha argument is %1%, but must be >= 0 !", alpha); + "Alpha argument is %1%, but must be > 0 !", alpha); return false; } return true; @@ -59,11 +59,11 @@ namespace boost template inline bool check_beta(const char* function, const RealType& beta, RealType* result) { - if(!(boost::math::isfinite)(beta) || (beta < 0)) + if(!(boost::math::isfinite)(beta) || (beta <= 0)) { *result = tools::domain_error( function, - "Beta argument is %1%, but must be >= 0 !", beta); + "Beta argument is %1%, but must be > 0 !", beta); return false; } return true; @@ -119,11 +119,11 @@ namespace boost template inline bool check_mean(const char* function, const RealType& mean, RealType* result) { - if(!(boost::math::isfinite)(mean) || (mean < 0)) + if(!(boost::math::isfinite)(mean) || (mean <= 0)) { *result = tools::domain_error( function, - "mean argument is %1%, but must be >= 0 !", mean); + "mean argument is %1%, but must be > 0 !", mean); return false; } return true; @@ -131,11 +131,11 @@ namespace boost template inline bool check_variance(const char* function, const RealType& variance, RealType* result) { - if(!(boost::math::isfinite)(variance) || (variance < 0)) + if(!(boost::math::isfinite)(variance) || (variance <= 0)) { *result = tools::domain_error( function, - "variance argument is %1%, but must be >= 0 !", variance); + "variance argument is %1%, but must be > 0 !", variance); return false; } return true; @@ -314,24 +314,23 @@ namespace boost return (2 * (b-a) * sqrt(a + b + 1)) / ((a + b + 2) * sqrt(a * b)); } // skewness - template - inline RealType kurtosis(const beta_distribution& dist) - { - RealType a = dist.alpha(); - RealType b = dist.beta(); - return 3 + 6 * a * a * a - a * a * (2 * b -1) + b * b * (b + 1) - 2 * a * b * (b + 2) / - a * b * (a + b + 2) * (a + b + 3); - } // kurtosis - template inline RealType kurtosis_excess(const beta_distribution& dist) { RealType a = dist.alpha(); RealType b = dist.beta(); - return 6 * a * a * a - a * a * (2 * b -1) + b * b * (b + 1) - 2 * a * b * (b + 2) / - a * b * (a + b + 2) * (a + b + 3); + RealType a_2 = a * a; + RealType n = 6 * (a_2 * a - a_2 * (2 * b - 1) + b * b * (b + 1) - 2 * a * b * (b + 2)); + RealType d = a * b * (a + b + 2) * (a + b + 3); + return n / d; } // kurtosis_excess + template + inline RealType kurtosis(const beta_distribution& dist) + { + return 3 + kurtosis_excess(dist); + } // kurtosis + template RealType pdf(const beta_distribution& dist, const RealType x) { // Probability Density/Mass Function. @@ -353,7 +352,7 @@ namespace boost return result; } using boost::math::beta; - return pow(x, (a-1)) * pow((1 - x), (b-1))/ beta(a, b); + return ibeta_derivative(a, b, x); } // pdf template @@ -418,7 +417,7 @@ namespace boost // Calculate cdf beta using the incomplete beta function. // Use of ibeta here prevents cancellation errors in calculating // 1 - x if x is very small, perhaps smaller than machine epsilon. - return 1 - ibeta(a, b, x); + return ibetac(a, b, x); } // beta cdf template diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index 15552ae99..1153de280 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -808,6 +808,15 @@ T ibeta_imp(T a, T b, T x, const L& l, bool inv, bool normalised) T fract; T y = 1 - x; + if(normalised) + { + // extend to a few very special cases: + if((a == 0) && (b != 0)) + return inv ? 0 : 1; + else if(b == 0) + return inv ? 1 : 0; + } + if(a <= 0) tools::domain_error(BOOST_CURRENT_FUNCTION, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a); if(b <= 0) @@ -1061,12 +1070,12 @@ T ibeta_derivative_imp(T a, T b, T x, const L& l) if(x == 0) { return (a > 1) ? 0 : - (a == 1) ? 1 : tools::overflow_error(BOOST_CURRENT_FUNCTION, 0); + (a == 1) ? 1 / boost::math::beta(a, b) : tools::overflow_error(BOOST_CURRENT_FUNCTION, 0); } else if(x == 1) { return (b > 1) ? 0 : - (b == 1) ? 1 : tools::overflow_error(BOOST_CURRENT_FUNCTION, 0); + (b == 1) ? 1 / boost::math::beta(a, b) : tools::overflow_error(BOOST_CURRENT_FUNCTION, 0); } // // Now the regular cases: diff --git a/test/test_beta_dist.cpp b/test/test_beta_dist.cpp index 7f37be653..78aa7369e 100644 --- a/test/test_beta_dist.cpp +++ b/test/test_beta_dist.cpp @@ -196,10 +196,11 @@ void test_spots(RealType) pdf(beta_distribution(static_cast(1), static_cast(1)), static_cast(0)), // x static_cast(1)); - BOOST_CHECK_EQUAL( + BOOST_CHECK_CLOSE_FRACTION( pdf(beta_distribution(static_cast(1), static_cast(1)), static_cast(0.5)), // x - static_cast(1)); + static_cast(1), + tolerance); BOOST_CHECK_EQUAL( beta_distribution(static_cast(1), static_cast(1)).alpha(), @@ -225,14 +226,14 @@ void test_spots(RealType) BOOST_CHECK_CLOSE_FRACTION( cdf(beta_distribution(static_cast(2), static_cast(2)), static_cast(0.1)), // x - static_cast(0.02800000000000000000000000000000000000000), // Seems exact. + static_cast(0.02800000000000000000000000000000000000000L), // Seems exact. // http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=BetaRegularized&ptype=0&z=0.1&a=2&b=2&digits=40 tolerance); BOOST_CHECK_CLOSE_FRACTION( cdf(beta_distribution(static_cast(2), static_cast(2)), static_cast(0.0001)), // x - static_cast(2.999800000000000000000000000000000000000e-8), + static_cast(2.999800000000000000000000000000000000000e-8L), // http://members.aol.com/iandjmsmith/BETAEX.HTM 2.9998000000004 // http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=BetaRegularized&ptype=0&z=0.0001&a=2&b=2&digits=40 tolerance); @@ -241,41 +242,42 @@ void test_spots(RealType) BOOST_CHECK_CLOSE_FRACTION( pdf(beta_distribution(static_cast(2), static_cast(2)), static_cast(0.0001)), // x - static_cast(0.0005999400000000004), // http://members.aol.com/iandjmsmith/BETAEX.HTM - tolerance); + static_cast(0.0005999400000000004L), // http://members.aol.com/iandjmsmith/BETAEX.HTM + // Slightly higher tolerance for real concept: + (std::numeric_limits::is_specialized ? 1 : 10) * tolerance); BOOST_CHECK_CLOSE_FRACTION( cdf(beta_distribution(static_cast(2), static_cast(2)), static_cast(0.9999)), // x - static_cast(0.999999970002), // http://members.aol.com/iandjmsmith/BETAEX.HTM + static_cast(0.999999970002L), // http://members.aol.com/iandjmsmith/BETAEX.HTM // Wolfram 0.9999999700020000000000000000000000000000 tolerance); BOOST_CHECK_CLOSE_FRACTION( cdf(beta_distribution(static_cast(0.5), static_cast(2)), static_cast(0.9)), // x - static_cast(0.9961174629530394895796514664963063381217), + static_cast(0.9961174629530394895796514664963063381217L), // Wolfram tolerance); BOOST_CHECK_CLOSE_FRACTION( cdf(beta_distribution(static_cast(0.5), static_cast(0.5)), static_cast(0.1)), // x - static_cast(0.2048327646991334516491978475505189480977), + static_cast(0.2048327646991334516491978475505189480977L), // Wolfram tolerance); BOOST_CHECK_CLOSE_FRACTION( cdf(beta_distribution(static_cast(0.5), static_cast(0.5)), static_cast(0.9)), // x - static_cast(0.7951672353008665483508021524494810519023), + static_cast(0.7951672353008665483508021524494810519023L), // Wolfram tolerance); BOOST_CHECK_CLOSE_FRACTION( quantile(beta_distribution(static_cast(0.5), static_cast(0.5)), - static_cast(0.7951672353008665483508021524494810519023)), // x + static_cast(0.7951672353008665483508021524494810519023L)), // x static_cast(0.9), // Wolfram tolerance); @@ -283,13 +285,13 @@ void test_spots(RealType) BOOST_CHECK_CLOSE_FRACTION( cdf(beta_distribution(static_cast(0.5), static_cast(0.5)), static_cast(0.6)), // x - static_cast(0.5640942168489749316118742861695149357858), + static_cast(0.5640942168489749316118742861695149357858L), // Wolfram tolerance); BOOST_CHECK_CLOSE_FRACTION( quantile(beta_distribution(static_cast(0.5), static_cast(0.5)), - static_cast(0.5640942168489749316118742861695149357858)), // x + static_cast(0.5640942168489749316118742861695149357858L)), // x static_cast(0.6), // Wolfram tolerance); @@ -298,13 +300,13 @@ void test_spots(RealType) BOOST_CHECK_CLOSE_FRACTION( cdf(beta_distribution(static_cast(2), static_cast(0.5)), static_cast(0.6)), // x - static_cast(0.1778078083562213736802876784474931812329), + static_cast(0.1778078083562213736802876784474931812329L), // Wolfram tolerance); BOOST_CHECK_CLOSE_FRACTION( quantile(beta_distribution(static_cast(2), static_cast(0.5)), - static_cast(0.1778078083562213736802876784474931812329)), // x + static_cast(0.1778078083562213736802876784474931812329L)), // x static_cast(0.6), // Wolfram tolerance); // gives @@ -326,13 +328,13 @@ void test_spots(RealType) BOOST_CHECK_CLOSE_FRACTION( cdf(complement(beta_distribution(static_cast(0.5), static_cast(0.5)), static_cast(0.1))), // complement of x - static_cast(0.7951672353008665483508021524494810519023), + static_cast(0.7951672353008665483508021524494810519023L), // Wolfram tolerance); BOOST_CHECK_CLOSE_FRACTION( quantile(beta_distribution(static_cast(2), static_cast(2)), - static_cast(0.0280000000000000000000000000000000000)), // x + static_cast(0.0280000000000000000000000000000000000L)), // x static_cast(0.1), // Wolfram tolerance); @@ -341,14 +343,14 @@ void test_spots(RealType) BOOST_CHECK_CLOSE_FRACTION( cdf(complement(beta_distribution(static_cast(2), static_cast(2)), static_cast(0.1))), // x - static_cast(0.9720000000000000000000000000000000000000), // Exact. + static_cast(0.9720000000000000000000000000000000000000L), // Exact. // Wolfram tolerance); BOOST_CHECK_CLOSE_FRACTION( pdf(beta_distribution(static_cast(2), static_cast(2)), static_cast(0.9999)), // x - static_cast(0.0005999399999999344), // http://members.aol.com/iandjmsmith/BETAEX.HTM + static_cast(0.0005999399999999344L), // http://members.aol.com/iandjmsmith/BETAEX.HTM tolerance*10); // Note less accurate. //void test_spot( @@ -374,8 +376,8 @@ void test_spots(RealType) static_cast(2), // alpha a static_cast(2), // beta b static_cast(0.1), // Probability p - static_cast(0.0280000000000000000000000000000000000), // Probability of result (CDF of beta), P - static_cast(1 - 0.0280000000000000000000000000000000000), // Complement of CDF Q = 1 - P + static_cast(0.0280000000000000000000000000000000000L), // Probability of result (CDF of beta), P + static_cast(1 - 0.0280000000000000000000000000000000000L), // Complement of CDF Q = 1 - P tolerance); // Test tolerance. @@ -399,56 +401,56 @@ void test_spots(RealType) static_cast(2), // alpha a static_cast(2), // beta b static_cast(0.01), // Probability p - static_cast(0.0002980000000000000000000000000000000000000), // Probability of result (CDF of beta), P - static_cast(1-0.0002980000000000000000000000000000000000000), // Complement of CDF Q = 1 - P + static_cast(0.0002980000000000000000000000000000000000000L), // Probability of result (CDF of beta), P + static_cast(1-0.0002980000000000000000000000000000000000000L), // Complement of CDF Q = 1 - P tolerance); // Test tolerance. test_spot( static_cast(2), // alpha a static_cast(2), // beta b static_cast(0.001), // Probability p - static_cast(2.998000000000000000000000000000000000000E-6), // Probability of result (CDF of beta), P - static_cast(1-2.998000000000000000000000000000000000000E-6), // Complement of CDF Q = 1 - P + static_cast(2.998000000000000000000000000000000000000E-6L), // Probability of result (CDF of beta), P + static_cast(1-2.998000000000000000000000000000000000000E-6L), // Complement of CDF Q = 1 - P tolerance); // Test tolerance. test_spot( static_cast(2), // alpha a static_cast(2), // beta b static_cast(0.0001), // Probability p - static_cast(2.999800000000000000000000000000000000000E-8), // Probability of result (CDF of beta), P - static_cast(1-2.999800000000000000000000000000000000000E-8), // Complement of CDF Q = 1 - P + static_cast(2.999800000000000000000000000000000000000E-8L), // Probability of result (CDF of beta), P + static_cast(1-2.999800000000000000000000000000000000000E-8L), // Complement of CDF Q = 1 - P tolerance); // Test tolerance. test_spot( static_cast(2), // alpha a static_cast(2), // beta b static_cast(0.99), // Probability p - static_cast(0.9997020000000000000000000000000000000000), // Probability of result (CDF of beta), P - static_cast(1-0.9997020000000000000000000000000000000000), // Complement of CDF Q = 1 - P + static_cast(0.9997020000000000000000000000000000000000L), // Probability of result (CDF of beta), P + static_cast(1-0.9997020000000000000000000000000000000000L), // Complement of CDF Q = 1 - P tolerance); // Test tolerance. test_spot( static_cast(0.5), // alpha a static_cast(2), // beta b static_cast(0.5), // Probability p - static_cast(0.8838834764831844055010554526310612991060), // Probability of result (CDF of beta), P - static_cast(1-0.8838834764831844055010554526310612991060), // Complement of CDF Q = 1 - P + static_cast(0.8838834764831844055010554526310612991060L), // Probability of result (CDF of beta), P + static_cast(1-0.8838834764831844055010554526310612991060L), // Complement of CDF Q = 1 - P tolerance); // Test tolerance. test_spot( static_cast(0.5), // alpha a static_cast(3.), // beta b static_cast(0.7), // Probability p - static_cast(0.9903963064097119299191611355232156905687), // Probability of result (CDF of beta), P - static_cast(1-0.9903963064097119299191611355232156905687), // Complement of CDF Q = 1 - P + static_cast(0.9903963064097119299191611355232156905687L), // Probability of result (CDF of beta), P + static_cast(1-0.9903963064097119299191611355232156905687L), // Complement of CDF Q = 1 - P tolerance); // Test tolerance. test_spot( static_cast(0.5), // alpha a static_cast(3.), // beta b static_cast(0.1), // Probability p - static_cast(0.5545844446520295253493059553548880128511), // Probability of result (CDF of beta), P - static_cast(1-0.5545844446520295253493059553548880128511), // Complement of CDF Q = 1 - P + static_cast(0.5545844446520295253493059553548880128511L), // Probability of result (CDF of beta), P + static_cast(1-0.5545844446520295253493059553548880128511L), // Complement of CDF Q = 1 - P tolerance); // Test tolerance. } // template void test_spots(RealType) @@ -462,99 +464,96 @@ int test_main(int, char* []) #endif - // Check that can generate beta distribution using one convenience methods: - beta_distribution<> mybeta11(1., 1.); // Using default RealType double. - // but that - //boost::math::beta mybeta1(1., 1.); // Using typedef fails. - // error C2039: 'beta' : is not a member of 'boost::math' + // Check that can generate beta distribution using one convenience methods: + beta_distribution<> mybeta11(1., 1.); // Using default RealType double. + // but that + // boost::math::beta mybeta1(1., 1.); // Using typedef fails. + // error C2039: 'beta' : is not a member of 'boost::math' - // Basic sanity-check spot values. + // Basic sanity-check spot values. - // Some simple checks using double only. - BOOST_CHECK_EQUAL(mybeta11.alpha(), 1); // - BOOST_CHECK_EQUAL(mybeta11.beta(), 1); - BOOST_CHECK_EQUAL(mean(mybeta11), 0.5); // 1 / (1 + 1) = 1/2 exactly - BOOST_CHECK_THROW(mode(mybeta11), std::domain_error); - beta_distribution<> mybeta22(2., 2.); // pdf is dome shape. - BOOST_CHECK_EQUAL(mode(mybeta22), 0.5); // 2-1 / (2+2-2) = 1/2 exactly. - beta_distribution<> mybetaH2(0.5, 2.); // - beta_distribution<> mybetaH3(0.5, 3.); // + // Some simple checks using double only. + BOOST_CHECK_EQUAL(mybeta11.alpha(), 1); // + BOOST_CHECK_EQUAL(mybeta11.beta(), 1); + BOOST_CHECK_EQUAL(mean(mybeta11), 0.5); // 1 / (1 + 1) = 1/2 exactly + BOOST_CHECK_THROW(mode(mybeta11), std::domain_error); + beta_distribution<> mybeta22(2., 2.); // pdf is dome shape. + BOOST_CHECK_EQUAL(mode(mybeta22), 0.5); // 2-1 / (2+2-2) = 1/2 exactly. + beta_distribution<> mybetaH2(0.5, 2.); // + beta_distribution<> mybetaH3(0.5, 3.); // - // Check a few values using double. - BOOST_CHECK_EQUAL(pdf(mybeta11, 1), 1); // is uniform unity over 0 to 1, - BOOST_CHECK_EQUAL(pdf(mybeta11, 0), 1); // including zero and unity. - BOOST_CHECK_EQUAL(pdf(mybeta11, 0.5), 1); - BOOST_CHECK_EQUAL(pdf(mybeta11, 0.0001), 1); - BOOST_CHECK_EQUAL(pdf(mybeta11, 0.9999), 1); - BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta11, 0.1), 0.1, std::numeric_limits::epsilon()); - BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta11, 0.5), 0.5, std::numeric_limits::epsilon()); - BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta11, 0.9), 0.9, std::numeric_limits::epsilon()); - BOOST_CHECK_EQUAL(cdf(mybeta11, 1), 1.); // Exact unity expected. + // Check a few values using double. + BOOST_CHECK_EQUAL(pdf(mybeta11, 1), 1); // is uniform unity over 0 to 1, + BOOST_CHECK_EQUAL(pdf(mybeta11, 0), 1); // including zero and unity. + // Although these next three have an exact result, internally they're + // *not* treated as special cases, and may be out by a couple of eps: + BOOST_CHECK_CLOSE_FRACTION(pdf(mybeta11, 0.5), 1.0, 5*std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE_FRACTION(pdf(mybeta11, 0.0001), 1.0, 5*std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE_FRACTION(pdf(mybeta11, 0.9999), 1.0, 5*std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta11, 0.1), 0.1, std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta11, 0.5), 0.5, std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta11, 0.9), 0.9, std::numeric_limits::epsilon()); + BOOST_CHECK_EQUAL(cdf(mybeta11, 1), 1.); // Exact unity expected. - double tol = std::numeric_limits::epsilon() * 10; - BOOST_CHECK_EQUAL(pdf(mybeta22, 1), 0); // is dome shape. - BOOST_CHECK_EQUAL(pdf(mybeta22, 0), 0); - BOOST_CHECK_CLOSE_FRACTION(pdf(mybeta22, 0.5), 1.5, tol); // top of dome, expect exactly 3/2. - BOOST_CHECK_CLOSE_FRACTION(pdf(mybeta22, 0.0001), 5.9994000000000E-4, tol); - BOOST_CHECK_CLOSE_FRACTION(pdf(mybeta22, 0.9999), 5.9994000000000E-4, tol*50); + double tol = std::numeric_limits::epsilon() * 10; + BOOST_CHECK_EQUAL(pdf(mybeta22, 1), 0); // is dome shape. + BOOST_CHECK_EQUAL(pdf(mybeta22, 0), 0); + BOOST_CHECK_CLOSE_FRACTION(pdf(mybeta22, 0.5), 1.5, tol); // top of dome, expect exactly 3/2. + BOOST_CHECK_CLOSE_FRACTION(pdf(mybeta22, 0.0001), 5.9994000000000E-4, tol); + BOOST_CHECK_CLOSE_FRACTION(pdf(mybeta22, 0.9999), 5.9994000000000E-4, tol*50); - BOOST_CHECK_EQUAL(cdf(mybeta22, 0.), 0); // cdf is a curved line from 0 to 1. - BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.1), 0.028000000000000, tol); - BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.5), 0.5, tol); - BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.9), 0.972000000000000, tol); - BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.0001), 2.999800000000000000000000000000000000000E-8, tol); - BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.001), 2.998000000000000000000000000000000000000E-6, tol); - BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.01), 0.0002980000000000000000000000000000000000000, tol); - BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.1), 0.02800000000000000000000000000000000000000, tol); // exact - BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.99), 0.9997020000000000000000000000000000000000, tol); + BOOST_CHECK_EQUAL(cdf(mybeta22, 0.), 0); // cdf is a curved line from 0 to 1. + BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.1), 0.028000000000000, tol); + BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.5), 0.5, tol); + BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.9), 0.972000000000000, tol); + BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.0001), 2.999800000000000000000000000000000000000E-8, tol); + BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.001), 2.998000000000000000000000000000000000000E-6, tol); + BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.01), 0.0002980000000000000000000000000000000000000, tol); + BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.1), 0.02800000000000000000000000000000000000000, tol); // exact + BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.99), 0.9997020000000000000000000000000000000000, tol); - BOOST_CHECK_EQUAL(cdf(mybeta22, 1), 1.); // Exact unity expected. + BOOST_CHECK_EQUAL(cdf(mybeta22, 1), 1.); // Exact unity expected. - // Complement + // Complement - BOOST_CHECK_CLOSE_FRACTION(cdf(complement(mybeta22, 0.9)), 0.028000000000000, tol); + BOOST_CHECK_CLOSE_FRACTION(cdf(complement(mybeta22, 0.9)), 0.028000000000000, tol); - // quantile. - BOOST_CHECK_CLOSE_FRACTION(quantile(mybeta22, 0.028), 0.1, tol); - BOOST_CHECK_CLOSE_FRACTION(quantile(complement(mybeta22, 1 - 0.028)), 0.1, tol); - BOOST_CHECK_EQUAL(kurtosis(mybeta11), 3+ kurtosis_excess(mybeta11)); // Check kurtosis_excess = kurtosis - 3; - BOOST_CHECK_CLOSE_FRACTION(variance(mybeta22), 0.05, tol); - BOOST_CHECK_CLOSE_FRACTION(mode(mybeta22), 0.5, tol); - BOOST_CHECK_CLOSE_FRACTION(mean(mybeta22), 0.5, tol); + // quantile. + BOOST_CHECK_CLOSE_FRACTION(quantile(mybeta22, 0.028), 0.1, tol); + BOOST_CHECK_CLOSE_FRACTION(quantile(complement(mybeta22, 1 - 0.028)), 0.1, tol); + BOOST_CHECK_EQUAL(kurtosis(mybeta11), 3+ kurtosis_excess(mybeta11)); // Check kurtosis_excess = kurtosis - 3; + BOOST_CHECK_CLOSE_FRACTION(variance(mybeta22), 0.05, tol); + BOOST_CHECK_CLOSE_FRACTION(mode(mybeta22), 0.5, tol); + BOOST_CHECK_CLOSE_FRACTION(mean(mybeta22), 0.5, tol); + BOOST_CHECK_CLOSE_FRACTION(skewness(mybeta22), 0.0, tol); + BOOST_CHECK_CLOSE_FRACTION(kurtosis_excess(mybeta22), -144.0 / 168, tol); + BOOST_CHECK_CLOSE_FRACTION(skewness(beta_distribution<>(3, 5)), 0.30983866769659335081434123198259, tol); + + BOOST_CHECK_EQUAL(beta_distribution::estimate_alpha(mean(mybeta22), variance(mybeta22)), mybeta22.alpha()); // mean, variance, probability. + BOOST_CHECK_EQUAL(beta_distribution::estimate_beta(mean(mybeta22), variance(mybeta22)), mybeta22.beta());// mean, variance, probability. + + BOOST_CHECK_CLOSE_FRACTION(mybeta22.estimate_alpha(mybeta22.beta(), 0.8, cdf(mybeta22, 0.8)), mybeta22.alpha(), tol); + BOOST_CHECK_CLOSE_FRACTION(mybeta22.estimate_beta(mybeta22.alpha(), 0.8, cdf(mybeta22, 0.8)), mybeta22.beta(), tol); - //cout << beta_distribution::estimate_alpha(mean(mybeta22), variance(mybeta22)) << endl; // 2 - //cout << beta_distribution::estimate_beta(mean(mybeta22), variance(mybeta22)) << endl; // 2 - BOOST_CHECK_EQUAL(beta_distribution::estimate_alpha(mean(mybeta22), variance(mybeta22)), mybeta22.alpha()); // mean, variance, probability. - BOOST_CHECK_EQUAL(beta_distribution::estimate_beta(mean(mybeta22), variance(mybeta22)), mybeta22.beta());// mean, variance, probability. - - // BOOST_CHECK_CLOSE_FRACTION(ibeta_inva(mybeta22.beta(), 0.8, cdf(mybeta22, 0.8)), mybeta22.alpha(), tol); - //BOOST_CHECK_CLOSE_FRACTION(ibeta_inva(mybeta22.beta(), 0.8, cdf(mybeta22, 0.8)), mybeta22.alpha(), tol); - BOOST_CHECK_CLOSE_FRACTION(mybeta22.estimate_alpha(mybeta22.beta(), 0.8, cdf(mybeta22, 0.8)), mybeta22.alpha(), tol); - BOOST_CHECK_CLOSE_FRACTION(mybeta22.estimate_beta(mybeta22.alpha(), 0.8, cdf(mybeta22, 0.8)), mybeta22.beta(), tol); + beta_distribution rcbeta22(2, 2); // Using RealType real_concept. + cout << "numeric_limits::is_specialized " << numeric_limits::is_specialized << endl; + cout << "numeric_limits::digits " << numeric_limits::digits << endl; + cout << "numeric_limits::digits10 " << numeric_limits::digits10 << endl; + cout << "numeric_limits::epsilon " << numeric_limits::epsilon() << endl; - beta_distribution rcbeta22(2, 2); // Using RealType real_concept. - cout << "numeric_limits::is_specialized " << numeric_limits::is_specialized << endl; - cout << "numeric_limits::digits " << numeric_limits::digits << endl; - cout << "numeric_limits::digits10 " << numeric_limits::digits10 << endl; - cout << "numeric_limits::epsilon " << numeric_limits::epsilon() << endl; - - // Tests for improvements to Boost.test display of errors. - //BOOST_CHECK_CLOSE_FRACTION(rcbeta22.alpha(), static_cast(2 + 3 * std::numeric_limits::epsilon()), 0); - //BOOST_CHECK_CLOSE_FRACTION(cdf(rcbeta22, 0.1), static_cast(0.028000000000009), 0); - - // (Parameter value, arbitrarily zero, only communicates the floating point type). - test_spots(0.0F); // Test float. - test_spots(0.0); // Test double. + // (Parameter value, arbitrarily zero, only communicates the floating point type). + test_spots(0.0F); // Test float. + test_spots(0.0); // Test double. #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS - test_spots(0.0L); // Test long double. + test_spots(0.0L); // Test long double. #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) - test_spots(boost::math::concepts::real_concept(0.)); // Test real concept. + test_spots(boost::math::concepts::real_concept(0.)); // Test real concept. #endif #endif - return 0; + return 0; } // int test_main(int, char* []) /*