# Statistics.cpp

Go to the documentation of this file.
```00001 /*
00002  * This program is free software; you can redistribute it and/or modify
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 }
```

SHOGUN Machine Learning Toolbox - Documentation