SHOGUN
v1.1.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Written (W) 2011 Heiko Strathmann 00008 * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society 00009 * 00010 * ALGLIB Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier under GPL2+ 00011 * http://www.alglib.net/ 00012 * See header file for which functions are taken from ALGLIB (with adjustments 00013 * for shogun) 00014 */ 00015 00016 #include <shogun/mathematics/Statistics.h> 00017 #include <shogun/mathematics/Math.h> 00018 00019 using namespace shogun; 00020 00021 float64_t CStatistics::mean(SGVector<float64_t> values) 00022 { 00023 ASSERT(values.vlen>0); 00024 ASSERT(values.vector); 00025 00026 float64_t sum=0; 00027 for (index_t i=0; i<values.vlen; ++i) 00028 sum+=values.vector[i]; 00029 00030 return sum/values.vlen; 00031 } 00032 00033 float64_t CStatistics::variance(SGVector<float64_t> values) 00034 { 00035 ASSERT(values.vlen>1); 00036 ASSERT(values.vector); 00037 00038 float64_t mean=CStatistics::mean(values); 00039 00040 float64_t sum_squared_diff=0; 00041 for (index_t i=0; i<values.vlen; ++i) 00042 sum_squared_diff+=CMath::pow(values.vector[i]-mean, 2); 00043 00044 return sum_squared_diff/(values.vlen-1); 00045 } 00046 00047 float64_t CStatistics::std_deviation(SGVector<float64_t> values) 00048 { 00049 return CMath::sqrt(variance(values)); 00050 } 00051 00052 float64_t CStatistics::confidence_intervals_mean(SGVector<float64_t> values, 00053 float64_t alpha, float64_t& conf_int_low, float64_t& conf_int_up) 00054 { 00055 ASSERT(values.vlen>1); 00056 ASSERT(values.vector); 00057 00058 /* using one sided student t distribution evaluation */ 00059 alpha=alpha/2; 00060 00061 /* degrees of freedom */ 00062 int32_t deg=values.vlen-1; 00063 00064 /* compute absolute value of t-value */ 00065 float64_t t=CMath::abs(inverse_student_t_distribution(deg, alpha)); 00066 00067 /* values for calculating confidence interval */ 00068 float64_t std_dev=std_deviation(values); 00069 float64_t mean=CStatistics::mean(values); 00070 00071 /* compute confidence interval */ 00072 float64_t interval=t*std_dev/CMath::sqrt((float64_t)values.vlen); 00073 conf_int_low=mean-interval; 00074 conf_int_up=mean+interval; 00075 00076 return mean; 00077 } 00078 00079 float64_t CStatistics::student_t_distribution(int32_t k, float64_t t) 00080 { 00081 float64_t x; 00082 float64_t rk; 00083 float64_t z; 00084 float64_t f; 00085 float64_t tz; 00086 float64_t p; 00087 float64_t xsqk; 00088 int32_t j; 00089 float64_t result; 00090 00091 ASSERT(k>0); 00092 if (t==0) 00093 { 00094 result=0.5; 00095 return result; 00096 } 00097 if (t<-2.0) 00098 { 00099 rk=k; 00100 z=rk/(rk+t*t); 00101 result=0.5*incomplete_beta(0.5*rk, 0.5, z); 00102 return result; 00103 } 00104 if (t<0) 00105 { 00106 x=-t; 00107 } 00108 else 00109 { 00110 x=t; 00111 } 00112 rk=k; 00113 z=1.0+x*x/rk; 00114 if (k%2!=0) 00115 { 00116 xsqk=x/CMath::sqrt(rk); 00117 p=CMath::atan(xsqk); 00118 if (k>1) 00119 { 00120 f=1.0; 00121 tz=1.0; 00122 j=3; 00123 while (j<=k-2&&tz/f>CMath::MACHINE_EPSILON) 00124 { 00125 tz=tz*((j-1)/(z*j)); 00126 f=f+tz; 00127 j=j+2; 00128 } 00129 p=p+f*xsqk/z; 00130 } 00131 p=p*2.0/CMath::PI; 00132 } 00133 else 00134 { 00135 f=1.0; 00136 tz=1.0; 00137 j=2; 00138 while (j<=k-2&&tz/f>CMath::MACHINE_EPSILON) 00139 { 00140 tz=tz*((j-1)/(z*j)); 00141 f=f+tz; 00142 j=j+2; 00143 } 00144 p=f*x/CMath::sqrt(z*rk); 00145 } 00146 if (t<0) 00147 { 00148 p=-p; 00149 } 00150 result=0.5+0.5*p; 00151 return result; 00152 } 00153 00154 float64_t CStatistics::incomplete_beta(float64_t a, float64_t b, float64_t x) 00155 { 00156 float64_t t; 00157 float64_t xc; 00158 float64_t w; 00159 float64_t y; 00160 int32_t flag; 00161 float64_t big; 00162 float64_t biginv; 00163 float64_t maxgam; 00164 float64_t minlog; 00165 float64_t maxlog; 00166 float64_t result; 00167 00168 big=4.503599627370496e15; 00169 biginv=2.22044604925031308085e-16; 00170 maxgam=171.624376956302725; 00171 minlog=CMath::log(CMath::MIN_REAL_NUMBER); 00172 maxlog=CMath::log(CMath::MAX_REAL_NUMBER); 00173 ASSERT(a>0&&b>0); 00174 ASSERT(x>=0&&x<=1); 00175 if (x==0) 00176 { 00177 result=0; 00178 return result; 00179 } 00180 if (x==1) 00181 { 00182 result=1; 00183 return result; 00184 } 00185 flag=0; 00186 if (b*x<=1.0&&x<=0.95) 00187 { 00188 result=ibetaf_incomplete_beta_ps(a, b, x, maxgam); 00189 return result; 00190 } 00191 w=1.0-x; 00192 if (x>a/(a+b)) 00193 { 00194 flag=1; 00195 t=a; 00196 a=b; 00197 b=t; 00198 xc=x; 00199 x=w; 00200 } 00201 else 00202 { 00203 xc=w; 00204 } 00205 if ((flag==1&&b*x<=1.0)&&x<=0.95) 00206 { 00207 t=ibetaf_incomplete_beta_ps(a, b, x, maxgam); 00208 if (t<=CMath::MACHINE_EPSILON) 00209 { 00210 result=1.0-CMath::MACHINE_EPSILON; 00211 } 00212 else 00213 { 00214 result=1.0-t; 00215 } 00216 return result; 00217 } 00218 y=x*(a+b-2.0)-(a-1.0); 00219 if (y<0.0) 00220 { 00221 w=ibetaf_incomplete_beta_fe(a, b, x, big, biginv); 00222 } 00223 else 00224 { 00225 w=ibetaf_incomplete_beta_fe2(a, b, x, big, biginv)/xc; 00226 } 00227 y=a*CMath::log(x); 00228 t=b*CMath::log(xc); 00229 if ((a+b<maxgam&&CMath::abs(y)<maxlog)&&CMath::abs(t)<maxlog) 00230 { 00231 t=CMath::pow(xc, b); 00232 t=t*CMath::pow(x, a); 00233 t=t/a; 00234 t=t*w; 00235 t=t*(CMath::tgamma(a+b)/(CMath::tgamma(a)*CMath::tgamma(b))); 00236 if (flag==1) 00237 { 00238 if (t<=CMath::MACHINE_EPSILON) 00239 { 00240 result=1.0-CMath::MACHINE_EPSILON; 00241 } 00242 else 00243 { 00244 result=1.0-t; 00245 } 00246 } 00247 else 00248 { 00249 result=t; 00250 } 00251 return result; 00252 } 00253 y=y+t+CMath::lgamma(a+b)-CMath::lgamma(a)-CMath::lgamma(b); 00254 y=y+CMath::log(w/a); 00255 if (y<minlog) 00256 { 00257 t=0.0; 00258 } 00259 else 00260 { 00261 t=CMath::exp(y); 00262 } 00263 if (flag==1) 00264 { 00265 if (t<=CMath::MACHINE_EPSILON) 00266 { 00267 t=1.0-CMath::MACHINE_EPSILON; 00268 } 00269 else 00270 { 00271 t=1.0-t; 00272 } 00273 } 00274 result=t; 00275 return result; 00276 } 00277 00278 float64_t CStatistics::ibetaf_incomplete_beta_ps(float64_t a, float64_t b, 00279 float64_t x, float64_t maxgam) 00280 { 00281 float64_t s; 00282 float64_t t; 00283 float64_t u; 00284 float64_t v; 00285 float64_t n; 00286 float64_t t1; 00287 float64_t z; 00288 float64_t ai; 00289 float64_t result; 00290 00291 ai=1.0/a; 00292 u=(1.0-b)*x; 00293 v=u/(a+1.0); 00294 t1=v; 00295 t=u; 00296 n=2.0; 00297 s=0.0; 00298 z=CMath::MACHINE_EPSILON*ai; 00299 while (CMath::abs(v)>z) 00300 { 00301 u=(n-b)*x/n; 00302 t=t*u; 00303 v=t/(a+n); 00304 s=s+v; 00305 n=n+1.0; 00306 } 00307 s=s+t1; 00308 s=s+ai; 00309 u=a*CMath::log(x); 00310 if (a+b<maxgam&&CMath::abs(u)<CMath::log(CMath::MAX_REAL_NUMBER)) 00311 { 00312 t=CMath::tgamma(a+b)/(CMath::tgamma(a)*CMath::tgamma(b)); 00313 s=s*t*CMath::pow(x, a); 00314 } 00315 else 00316 { 00317 t=CMath::lgamma(a+b)-CMath::lgamma(a)-CMath::lgamma(b)+u+CMath::log(s); 00318 if (t<CMath::log(CMath::MIN_REAL_NUMBER)) 00319 { 00320 s=0.0; 00321 } 00322 else 00323 { 00324 s=CMath::exp(t); 00325 } 00326 } 00327 result=s; 00328 return result; 00329 } 00330 00331 float64_t CStatistics::ibetaf_incomplete_beta_fe2(float64_t a, float64_t b, 00332 float64_t x, float64_t big, float64_t biginv) 00333 { 00334 float64_t xk; 00335 float64_t pk; 00336 float64_t pkm1; 00337 float64_t pkm2; 00338 float64_t qk; 00339 float64_t qkm1; 00340 float64_t qkm2; 00341 float64_t k1; 00342 float64_t k2; 00343 float64_t k3; 00344 float64_t k4; 00345 float64_t k5; 00346 float64_t k6; 00347 float64_t k7; 00348 float64_t k8; 00349 float64_t r; 00350 float64_t t; 00351 float64_t ans; 00352 float64_t z; 00353 float64_t thresh; 00354 int32_t n; 00355 float64_t result; 00356 00357 k1=a; 00358 k2=b-1.0; 00359 k3=a; 00360 k4=a+1.0; 00361 k5=1.0; 00362 k6=a+b; 00363 k7=a+1.0; 00364 k8=a+2.0; 00365 pkm2=0.0; 00366 qkm2=1.0; 00367 pkm1=1.0; 00368 qkm1=1.0; 00369 z=x/(1.0-x); 00370 ans=1.0; 00371 r=1.0; 00372 n=0; 00373 thresh=3.0*CMath::MACHINE_EPSILON; 00374 do 00375 { 00376 xk=-z*k1*k2/(k3*k4); 00377 pk=pkm1+pkm2*xk; 00378 qk=qkm1+qkm2*xk; 00379 pkm2=pkm1; 00380 pkm1=pk; 00381 qkm2=qkm1; 00382 qkm1=qk; 00383 xk=z*k5*k6/(k7*k8); 00384 pk=pkm1+pkm2*xk; 00385 qk=qkm1+qkm2*xk; 00386 pkm2=pkm1; 00387 pkm1=pk; 00388 qkm2=qkm1; 00389 qkm1=qk; 00390 if (qk!=0) 00391 { 00392 r=pk/qk; 00393 } 00394 if (r!=0) 00395 { 00396 t=CMath::abs((ans-r)/r); 00397 ans=r; 00398 } 00399 else 00400 { 00401 t=1.0; 00402 } 00403 if (t<thresh) 00404 { 00405 break; 00406 } 00407 k1=k1+1.0; 00408 k2=k2-1.0; 00409 k3=k3+2.0; 00410 k4=k4+2.0; 00411 k5=k5+1.0; 00412 k6=k6+1.0; 00413 k7=k7+2.0; 00414 k8=k8+2.0; 00415 if (CMath::abs(qk)+CMath::abs(pk)>big) 00416 { 00417 pkm2=pkm2*biginv; 00418 pkm1=pkm1*biginv; 00419 qkm2=qkm2*biginv; 00420 qkm1=qkm1*biginv; 00421 } 00422 if (CMath::abs(qk)<biginv||CMath::abs(pk)<biginv) 00423 { 00424 pkm2=pkm2*big; 00425 pkm1=pkm1*big; 00426 qkm2=qkm2*big; 00427 qkm1=qkm1*big; 00428 } 00429 n=n+1; 00430 } while (n!=300); 00431 result=ans; 00432 return result; 00433 } 00434 00435 float64_t CStatistics::ibetaf_incomplete_beta_fe(float64_t a, float64_t b, 00436 float64_t x, float64_t big, float64_t biginv) 00437 { 00438 float64_t xk; 00439 float64_t pk; 00440 float64_t pkm1; 00441 float64_t pkm2; 00442 float64_t qk; 00443 float64_t qkm1; 00444 float64_t qkm2; 00445 float64_t k1; 00446 float64_t k2; 00447 float64_t k3; 00448 float64_t k4; 00449 float64_t k5; 00450 float64_t k6; 00451 float64_t k7; 00452 float64_t k8; 00453 float64_t r; 00454 float64_t t; 00455 float64_t ans; 00456 float64_t thresh; 00457 int32_t n; 00458 float64_t result; 00459 00460 k1=a; 00461 k2=a+b; 00462 k3=a; 00463 k4=a+1.0; 00464 k5=1.0; 00465 k6=b-1.0; 00466 k7=k4; 00467 k8=a+2.0; 00468 pkm2=0.0; 00469 qkm2=1.0; 00470 pkm1=1.0; 00471 qkm1=1.0; 00472 ans=1.0; 00473 r=1.0; 00474 n=0; 00475 thresh=3.0*CMath::MACHINE_EPSILON; 00476 do 00477 { 00478 xk=-x*k1*k2/(k3*k4); 00479 pk=pkm1+pkm2*xk; 00480 qk=qkm1+qkm2*xk; 00481 pkm2=pkm1; 00482 pkm1=pk; 00483 qkm2=qkm1; 00484 qkm1=qk; 00485 xk=x*k5*k6/(k7*k8); 00486 pk=pkm1+pkm2*xk; 00487 qk=qkm1+qkm2*xk; 00488 pkm2=pkm1; 00489 pkm1=pk; 00490 qkm2=qkm1; 00491 qkm1=qk; 00492 if (qk!=0) 00493 { 00494 r=pk/qk; 00495 } 00496 if (r!=0) 00497 { 00498 t=CMath::abs((ans-r)/r); 00499 ans=r; 00500 } 00501 else 00502 { 00503 t=1.0; 00504 } 00505 if (t<thresh) 00506 { 00507 break; 00508 } 00509 k1=k1+1.0; 00510 k2=k2+1.0; 00511 k3=k3+2.0; 00512 k4=k4+2.0; 00513 k5=k5+1.0; 00514 k6=k6-1.0; 00515 k7=k7+2.0; 00516 k8=k8+2.0; 00517 if (CMath::abs(qk)+CMath::abs(pk)>big) 00518 { 00519 pkm2=pkm2*biginv; 00520 pkm1=pkm1*biginv; 00521 qkm2=qkm2*biginv; 00522 qkm1=qkm1*biginv; 00523 } 00524 if (CMath::abs(qk)<biginv||CMath::abs(pk)<biginv) 00525 { 00526 pkm2=pkm2*big; 00527 pkm1=pkm1*big; 00528 qkm2=qkm2*big; 00529 qkm1=qkm1*big; 00530 } 00531 n=n+1; 00532 } while (n!=300); 00533 result=ans; 00534 return result; 00535 } 00536 00537 float64_t CStatistics::inverse_student_t_distribution(int32_t k, float64_t p) 00538 { 00539 float64_t t; 00540 float64_t rk; 00541 float64_t z; 00542 int32_t rflg; 00543 float64_t result; 00544 00545 ASSERT(k>0&&p>0&&p<1); 00546 rk=k; 00547 if (p>0.25&&p<0.75) 00548 { 00549 if (p==0.5) 00550 { 00551 result=0; 00552 return result; 00553 } 00554 z=1.0-2.0*p; 00555 z=inverse_incomplete_beta(0.5, 0.5*rk, CMath::abs(z)); 00556 t=CMath::sqrt(rk*z/(1.0-z)); 00557 if (p<0.5) 00558 { 00559 t=-t; 00560 } 00561 result=t; 00562 return result; 00563 } 00564 rflg=-1; 00565 if (p>=0.5) 00566 { 00567 p=1.0-p; 00568 rflg=1; 00569 } 00570 z=inverse_incomplete_beta(0.5*rk, 0.5, 2.0*p); 00571 if (CMath::MAX_REAL_NUMBER*z<rk) 00572 { 00573 result=rflg*CMath::MAX_REAL_NUMBER; 00574 return result; 00575 } 00576 t=CMath::sqrt(rk/z-rk); 00577 result=rflg*t; 00578 return result; 00579 } 00580 00581 float64_t CStatistics::inverse_incomplete_beta(float64_t a, float64_t b, 00582 float64_t y) 00583 { 00584 float64_t aaa; 00585 float64_t bbb; 00586 float64_t y0; 00587 float64_t d; 00588 float64_t yyy; 00589 float64_t x; 00590 float64_t x0; 00591 float64_t x1; 00592 float64_t lgm; 00593 float64_t yp; 00594 float64_t di; 00595 float64_t dithresh; 00596 float64_t yl; 00597 float64_t yh; 00598 float64_t xt; 00599 int32_t i; 00600 int32_t rflg; 00601 int32_t dir; 00602 int32_t nflg; 00603 int32_t mainlooppos; 00604 int32_t ihalve; 00605 int32_t ihalvecycle; 00606 int32_t newt; 00607 int32_t newtcycle; 00608 int32_t breaknewtcycle; 00609 int32_t breakihalvecycle; 00610 float64_t result; 00611 00612 i=0; 00613 ASSERT(y>=0&&y<=1); 00614 00615 /* 00616 * special cases 00617 */ 00618 if (y==0) 00619 { 00620 result=0; 00621 return result; 00622 } 00623 if (y==1.0) 00624 { 00625 result=1; 00626 return result; 00627 } 00628 00629 /* 00630 * these initializations are not really necessary, 00631 * but without them compiler complains about 'possibly uninitialized variables'. 00632 */ 00633 dithresh=0; 00634 rflg=0; 00635 aaa=0; 00636 bbb=0; 00637 y0=0; 00638 x=0; 00639 yyy=0; 00640 lgm=0; 00641 dir=0; 00642 di=0; 00643 00644 /* 00645 * normal initializations 00646 */ 00647 x0=0.0; 00648 yl=0.0; 00649 x1=1.0; 00650 yh=1.0; 00651 nflg=0; 00652 mainlooppos=0; 00653 ihalve=1; 00654 ihalvecycle=2; 00655 newt=3; 00656 newtcycle=4; 00657 breaknewtcycle=5; 00658 breakihalvecycle=6; 00659 00660 /* 00661 * main loop 00662 */ 00663 for (;;) 00664 { 00665 00666 /* 00667 * start 00668 */ 00669 if (mainlooppos==0) 00670 { 00671 if (a<=1.0||b<=1.0) 00672 { 00673 dithresh=1.0e-6; 00674 rflg=0; 00675 aaa=a; 00676 bbb=b; 00677 y0=y; 00678 x=aaa/(aaa+bbb); 00679 yyy=incomplete_beta(aaa, bbb, x); 00680 mainlooppos=ihalve; 00681 continue; 00682 } 00683 else 00684 { 00685 dithresh=1.0e-4; 00686 } 00687 yp=-inverse_normal_distribution(y); 00688 if (y>0.5) 00689 { 00690 rflg=1; 00691 aaa=b; 00692 bbb=a; 00693 y0=1.0-y; 00694 yp=-yp; 00695 } 00696 else 00697 { 00698 rflg=0; 00699 aaa=a; 00700 bbb=b; 00701 y0=y; 00702 } 00703 lgm=(yp*yp-3.0)/6.0; 00704 x=2.0/(1.0/(2.0*aaa-1.0)+1.0/(2.0*bbb-1.0)); 00705 d=yp*CMath::sqrt(x+lgm)/x 00706 -(1.0/(2.0*bbb-1.0)-1.0/(2.0*aaa-1.0)) 00707 *(lgm+5.0/6.0-2.0/(3.0*x)); 00708 d=2.0*d; 00709 if (d<CMath::log(CMath::MIN_REAL_NUMBER)) 00710 { 00711 x=0; 00712 break; 00713 } 00714 x=aaa/(aaa+bbb*CMath::exp(d)); 00715 yyy=incomplete_beta(aaa, bbb, x); 00716 yp=(yyy-y0)/y0; 00717 if (CMath::abs(yp)<0.2) 00718 { 00719 mainlooppos=newt; 00720 continue; 00721 } 00722 mainlooppos=ihalve; 00723 continue; 00724 } 00725 00726 /* 00727 * ihalve 00728 */ 00729 if (mainlooppos==ihalve) 00730 { 00731 dir=0; 00732 di=0.5; 00733 i=0; 00734 mainlooppos=ihalvecycle; 00735 continue; 00736 } 00737 00738 /* 00739 * ihalvecycle 00740 */ 00741 if (mainlooppos==ihalvecycle) 00742 { 00743 if (i<=99) 00744 { 00745 if (i!=0) 00746 { 00747 x=x0+di*(x1-x0); 00748 if (x==1.0) 00749 { 00750 x=1.0-CMath::MACHINE_EPSILON; 00751 } 00752 if (x==0.0) 00753 { 00754 di=0.5; 00755 x=x0+di*(x1-x0); 00756 if (x==0.0) 00757 { 00758 break; 00759 } 00760 } 00761 yyy=incomplete_beta(aaa, bbb, x); 00762 yp=(x1-x0)/(x1+x0); 00763 if (CMath::abs(yp)<dithresh) 00764 { 00765 mainlooppos=newt; 00766 continue; 00767 } 00768 yp=(yyy-y0)/y0; 00769 if (CMath::abs(yp)<dithresh) 00770 { 00771 mainlooppos=newt; 00772 continue; 00773 } 00774 } 00775 if (yyy<y0) 00776 { 00777 x0=x; 00778 yl=yyy; 00779 if (dir<0) 00780 { 00781 dir=0; 00782 di=0.5; 00783 } 00784 else 00785 { 00786 if (dir>3) 00787 { 00788 di=1.0-(1.0-di)*(1.0-di); 00789 } 00790 else 00791 { 00792 if (dir>1) 00793 { 00794 di=0.5*di+0.5; 00795 } 00796 else 00797 { 00798 di=(y0-yyy)/(yh-yl); 00799 } 00800 } 00801 } 00802 dir=dir+1; 00803 if (x0>0.75) 00804 { 00805 if (rflg==1) 00806 { 00807 rflg=0; 00808 aaa=a; 00809 bbb=b; 00810 y0=y; 00811 } 00812 else 00813 { 00814 rflg=1; 00815 aaa=b; 00816 bbb=a; 00817 y0=1.0-y; 00818 } 00819 x=1.0-x; 00820 yyy=incomplete_beta(aaa, bbb, x); 00821 x0=0.0; 00822 yl=0.0; 00823 x1=1.0; 00824 yh=1.0; 00825 mainlooppos=ihalve; 00826 continue; 00827 } 00828 } 00829 else 00830 { 00831 x1=x; 00832 if (rflg==1&&x1<CMath::MACHINE_EPSILON) 00833 { 00834 x=0.0; 00835 break; 00836 } 00837 yh=yyy; 00838 if (dir>0) 00839 { 00840 dir=0; 00841 di=0.5; 00842 } 00843 else 00844 { 00845 if (dir<-3) 00846 { 00847 di=di*di; 00848 } 00849 else 00850 { 00851 if (dir<-1) 00852 { 00853 di=0.5*di; 00854 } 00855 else 00856 { 00857 di=(yyy-y0)/(yh-yl); 00858 } 00859 } 00860 } 00861 dir=dir-1; 00862 } 00863 i=i+1; 00864 mainlooppos=ihalvecycle; 00865 continue; 00866 } 00867 else 00868 { 00869 mainlooppos=breakihalvecycle; 00870 continue; 00871 } 00872 } 00873 00874 /* 00875 * breakihalvecycle 00876 */ 00877 if (mainlooppos==breakihalvecycle) 00878 { 00879 if (x0>=1.0) 00880 { 00881 x=1.0-CMath::MACHINE_EPSILON; 00882 break; 00883 } 00884 if (x<=0.0) 00885 { 00886 x=0.0; 00887 break; 00888 } 00889 mainlooppos=newt; 00890 continue; 00891 } 00892 00893 /* 00894 * newt 00895 */ 00896 if (mainlooppos==newt) 00897 { 00898 if (nflg!=0) 00899 { 00900 break; 00901 } 00902 nflg=1; 00903 lgm=CMath::lgamma(aaa+bbb)-CMath::lgamma(aaa)-CMath::lgamma(bbb); 00904 i=0; 00905 mainlooppos=newtcycle; 00906 continue; 00907 } 00908 00909 /* 00910 * newtcycle 00911 */ 00912 if (mainlooppos==newtcycle) 00913 { 00914 if (i<=7) 00915 { 00916 if (i!=0) 00917 { 00918 yyy=incomplete_beta(aaa, bbb, x); 00919 } 00920 if (yyy<yl) 00921 { 00922 x=x0; 00923 yyy=yl; 00924 } 00925 else 00926 { 00927 if (yyy>yh) 00928 { 00929 x=x1; 00930 yyy=yh; 00931 } 00932 else 00933 { 00934 if (yyy<y0) 00935 { 00936 x0=x; 00937 yl=yyy; 00938 } 00939 else 00940 { 00941 x1=x; 00942 yh=yyy; 00943 } 00944 } 00945 } 00946 if (x==1.0||x==0.0) 00947 { 00948 mainlooppos=breaknewtcycle; 00949 continue; 00950 } 00951 d=(aaa-1.0)*CMath::log(x)+(bbb-1.0)*CMath::log(1.0-x)+lgm; 00952 if (d<CMath::log(CMath::MIN_REAL_NUMBER)) 00953 { 00954 break; 00955 } 00956 if (d>CMath::log(CMath::MAX_REAL_NUMBER)) 00957 { 00958 mainlooppos=breaknewtcycle; 00959 continue; 00960 } 00961 d=CMath::exp(d); 00962 d=(yyy-y0)/d; 00963 xt=x-d; 00964 if (xt<=x0) 00965 { 00966 yyy=(x-x0)/(x1-x0); 00967 xt=x0+0.5*yyy*(x-x0); 00968 if (xt<=0.0) 00969 { 00970 mainlooppos=breaknewtcycle; 00971 continue; 00972 } 00973 } 00974 if (xt>=x1) 00975 { 00976 yyy=(x1-x)/(x1-x0); 00977 xt=x1-0.5*yyy*(x1-x); 00978 if (xt>=1.0) 00979 { 00980 mainlooppos=breaknewtcycle; 00981 continue; 00982 } 00983 } 00984 x=xt; 00985 if (CMath::abs(d/x)<128.0*CMath::MACHINE_EPSILON) 00986 { 00987 break; 00988 } 00989 i=i+1; 00990 mainlooppos=newtcycle; 00991 continue; 00992 } 00993 else 00994 { 00995 mainlooppos=breaknewtcycle; 00996 continue; 00997 } 00998 } 00999 01000 /* 01001 * breaknewtcycle 01002 */ 01003 if (mainlooppos==breaknewtcycle) 01004 { 01005 dithresh=256.0*CMath::MACHINE_EPSILON; 01006 mainlooppos=ihalve; 01007 continue; 01008 } 01009 } 01010 01011 /* 01012 * done 01013 */ 01014 if (rflg!=0) 01015 { 01016 if (x<=CMath::MACHINE_EPSILON) 01017 { 01018 x=1.0-CMath::MACHINE_EPSILON; 01019 } 01020 else 01021 { 01022 x=1.0-x; 01023 } 01024 } 01025 result=x; 01026 return result; 01027 } 01028 01029 float64_t CStatistics::inverse_normal_distribution(float64_t y0) 01030 { 01031 float64_t expm2; 01032 float64_t s2pi; 01033 float64_t x; 01034 float64_t y; 01035 float64_t z; 01036 float64_t y2; 01037 float64_t x0; 01038 float64_t x1; 01039 int32_t code; 01040 float64_t p0; 01041 float64_t q0; 01042 float64_t p1; 01043 float64_t q1; 01044 float64_t p2; 01045 float64_t q2; 01046 float64_t result; 01047 01048 expm2=0.13533528323661269189; 01049 s2pi=2.50662827463100050242; 01050 if (y0<=0) 01051 { 01052 result=-CMath::MAX_REAL_NUMBER; 01053 return result; 01054 } 01055 if (y0>=1) 01056 { 01057 result=CMath::MAX_REAL_NUMBER; 01058 return result; 01059 } 01060 code=1; 01061 y=y0; 01062 if (y>1.0-expm2) 01063 { 01064 y=1.0-y; 01065 code=0; 01066 } 01067 if (y>expm2) 01068 { 01069 y=y-0.5; 01070 y2=y*y; 01071 p0=-59.9633501014107895267; 01072 p0=98.0010754185999661536+y2*p0; 01073 p0=-56.6762857469070293439+y2*p0; 01074 p0=13.9312609387279679503+y2*p0; 01075 p0=-1.23916583867381258016+y2*p0; 01076 q0=1; 01077 q0=1.95448858338141759834+y2*q0; 01078 q0=4.67627912898881538453+y2*q0; 01079 q0=86.3602421390890590575+y2*q0; 01080 q0=-225.462687854119370527+y2*q0; 01081 q0=200.260212380060660359+y2*q0; 01082 q0=-82.0372256168333339912+y2*q0; 01083 q0=15.9056225126211695515+y2*q0; 01084 q0=-1.18331621121330003142+y2*q0; 01085 x=y+y*y2*p0/q0; 01086 x=x*s2pi; 01087 result=x; 01088 return result; 01089 } 01090 x=CMath::sqrt(-2.0*CMath::log(y)); 01091 x0=x-CMath::log(x)/x; 01092 z=1.0/x; 01093 if (x<8.0) 01094 { 01095 p1=4.05544892305962419923; 01096 p1=31.5251094599893866154+z*p1; 01097 p1=57.1628192246421288162+z*p1; 01098 p1=44.0805073893200834700+z*p1; 01099 p1=14.6849561928858024014+z*p1; 01100 p1=2.18663306850790267539+z*p1; 01101 p1=-1.40256079171354495875*0.1+z*p1; 01102 p1=-3.50424626827848203418*0.01+z*p1; 01103 p1=-8.57456785154685413611*0.0001+z*p1; 01104 q1=1; 01105 q1=15.7799883256466749731+z*q1; 01106 q1=45.3907635128879210584+z*q1; 01107 q1=41.3172038254672030440+z*q1; 01108 q1=15.0425385692907503408+z*q1; 01109 q1=2.50464946208309415979+z*q1; 01110 q1=-1.42182922854787788574*0.1+z*q1; 01111 q1=-3.80806407691578277194*0.01+z*q1; 01112 q1=-9.33259480895457427372*0.0001+z*q1; 01113 x1=z*p1/q1; 01114 } 01115 else 01116 { 01117 p2=3.23774891776946035970; 01118 p2=6.91522889068984211695+z*p2; 01119 p2=3.93881025292474443415+z*p2; 01120 p2=1.33303460815807542389+z*p2; 01121 p2=2.01485389549179081538*0.1+z*p2; 01122 p2=1.23716634817820021358*0.01+z*p2; 01123 p2=3.01581553508235416007*0.0001+z*p2; 01124 p2=2.65806974686737550832*0.000001+z*p2; 01125 p2=6.23974539184983293730*0.000000001+z*p2; 01126 q2=1; 01127 q2=6.02427039364742014255+z*q2; 01128 q2=3.67983563856160859403+z*q2; 01129 q2=1.37702099489081330271+z*q2; 01130 q2=2.16236993594496635890*0.1+z*q2; 01131 q2=1.34204006088543189037*0.01+z*q2; 01132 q2=3.28014464682127739104*0.0001+z*q2; 01133 q2=2.89247864745380683936*0.000001+z*q2; 01134 q2=6.79019408009981274425*0.000000001+z*q2; 01135 x1=z*p2/q2; 01136 } 01137 x=x0-x1; 01138 if (code!=0) 01139 { 01140 x=-x; 01141 } 01142 result=x; 01143 return result; 01144 }