PLearn 0.1
random.cc
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // PLearn (A C++ Machine Learning Library)
00004 // Copyright (C) 1998 Pascal Vincent
00005 // Copyright (C) 1999-2002 Pascal Vincent, Yoshua Bengio and University of Montreal
00006 //
00007 
00008 // Redistribution and use in source and binary forms, with or without
00009 // modification, are permitted provided that the following conditions are met:
00010 // 
00011 //  1. Redistributions of source code must retain the above copyright
00012 //     notice, this list of conditions and the following disclaimer.
00013 // 
00014 //  2. Redistributions in binary form must reproduce the above copyright
00015 //     notice, this list of conditions and the following disclaimer in the
00016 //     documentation and/or other materials provided with the distribution.
00017 // 
00018 //  3. The name of the authors may not be used to endorse or promote
00019 //     products derived from this software without specific prior written
00020 //     permission.
00021 // 
00022 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
00023 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
00024 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
00025 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00026 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
00027 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00028 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00029 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00030 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00031 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00032 // 
00033 // This file is part of the PLearn library. For more information on the PLearn
00034 // library, go to the PLearn Web site at www.plearn.org
00035 
00036 /* *******************************************************      
00037  * $Id: random.cc 8210 2007-10-24 19:12:31Z nouiz $
00038  ******************************************************* */
00039 
00040 extern "C" {
00041 #include <ctime>
00042 }
00043 
00044 #include <plearn/base/general.h>
00045 #include "random.h"
00046 
00047 namespace PLearn {
00048 using namespace std;
00049 
00050 
00051 /*  
00052     The static data to store the seed used by the random number generators.
00053 */
00054 
00055 static int32_t  the_seed=0;
00056 static int      iset=0;
00057 static real gset;
00058 
00059 /*  
00060     Special functions.
00061     =================
00062 */
00063 
00064 
00065 real  log_gamma(real xx)
00066 {
00067     double x,y,tmp,ser;
00068     static double gamma_coeffs[6]={ 76.18009172947146     ,
00069                                     -86.50532032941677     ,
00070                                     24.01409824083091     ,
00071                                     -1.231739572450155    ,
00072                                     0.1208650973866179e-2,
00073                                     -0.5395239384953e-5   };
00074     int j;
00075 
00076     y=x=xx;
00077     tmp=x+5.5;
00078     tmp -= (x+0.5)*pl_log(tmp);
00079     ser=1.000000000190015;
00080     for (j=0;j<=5;j++) ser += gamma_coeffs[j]/++y;
00081     return -tmp+pl_log(2.5066282746310005*ser/x);
00082 }
00083 
00084 real log_beta(real x, real y)
00085 {
00086     return log_gamma(x) + log_gamma(y) - log_gamma(x+y);
00087 }
00088 
00089 real incomplete_beta_continued_fraction(real z, real x, real y)
00090 {
00091     real x_minus_1 = x-1;
00092     real x_plus_1 = x+1;
00093     real x_plus_y = x+y;
00094     real denom = -z*x_plus_y/x_plus_1+1;
00095     if (fabs(denom)<1e-35) {
00096         denom=1e-35;
00097     }
00098     real rat1=1/denom;
00099     real rat2=1.0;
00100     real frac=rat1;
00101     for (int k=1;k<100;k++)
00102     {
00103         real f=z*k*(y-k)/((x+2*k)*(x_minus_1+2*k));
00104         rat2 = f/rat2 + 1;
00105         rat1 = rat1*f+1;
00106         if (fabs(rat1)<1e-35) {
00107             rat1=1e-35;
00108         }
00109         if (fabs(rat2)<1e-35) {
00110             rat2=1e-35;
00111         }
00112         rat1=1/rat1;
00113         frac *= rat1*rat2;
00114 
00115         f=-z*(x+k)*(x_plus_y+k)/((x_plus_1+2*k)*(x+2*k));
00116         rat2 = f/rat2+ 1;
00117         rat1 = rat1*f+1;
00118         if (fabs(rat1)<1e-35) {
00119             rat1=1e-35;
00120         }
00121         if (fabs(rat2)<1e-35) {
00122             rat2=1e-35;
00123         }
00124         rat1=1/rat1;
00125 
00126         real delta = rat1*rat2;
00127         frac *= delta;
00128         // stopping criterion
00129         if (fabs(1-delta) < 2e-7) {
00130             return frac;
00131         }
00132     }
00133     // If that happens, increase the number of k iterations or increase
00134     // the stopping criterion tolerance.
00135     PLWARNING("incomplete_beta_continued_fraction: insufficient precision!"); 
00136     return frac;
00137 }
00138 
00139 real incomplete_beta(real z, real x, real y)
00140 {
00141     if (z>1 || z<0) PLERROR("incomplete_beta(z,x,y): z =%f must be in [0,1]",z);
00142     real coeff = 0;
00143     if (z>0 && z<1) coeff = exp(x*pl_log(z)+y*pl_log(1.-z)-log_beta(x,y));
00144     if (z*(x+y+2)<x+1) {
00145         return coeff*incomplete_beta_continued_fraction(z,x,y)/x;
00146     }
00147     return 1-coeff*incomplete_beta_continued_fraction(1-z,y,x)/y;
00148 }
00149 
00150 real student_t_cdf(real t, int nb_degrees_of_freedom)
00151 {
00152     real p_t = 0.5*incomplete_beta(nb_degrees_of_freedom/(nb_degrees_of_freedom+t*t),0.5*nb_degrees_of_freedom,0.5);
00153     //real p_t = 0.5*incbet(0.5*nb_degrees_of_freedom,0.5,nb_degrees_of_freedom/(nb_degrees_of_freedom+t*t));
00154 #ifdef BOUNDCHECK
00155     if (p_t < 0) {
00156         PLERROR("Bug in incomplete_beta : returned a negative p_t !\n- p_t = %f\n- degrees of freedom = %d\n- t = %f",
00157                 p_t, nb_degrees_of_freedom, t);
00158     }
00159 #endif
00160     if (t>0)
00161         return 1.0 - p_t;
00162     else
00163         return p_t;
00164 }
00165 
00166 /*   
00167      Utilities for random numbers generation. 
00168      =======================================
00169 */
00170 
00171 /*  
00172     manual_seed(): gives a seed for random number generators.
00173 
00174     Rem: - The stored value is negative.
00175 */
00176 
00177 void  manual_seed(int32_t x)
00178 {
00179     the_seed = - labs(x);
00180     iset     = 0;
00181 }
00182 
00183 void  seed()
00184 {
00185     time_t  ltime;
00186     struct  tm *today;
00187     time(&ltime);
00188     today = localtime(&ltime);
00189     manual_seed((int32_t)today->tm_sec+
00190                 60*today->tm_min+
00191                 60*60*today->tm_hour+
00192                 60*60*24*today->tm_mday);
00193 }
00194 
00195 int32_t get_seed()
00196 {
00197     int32_t seed = the_seed;
00198     return seed;
00199 }
00200 
00201 /*  
00202     Constants used by the 'numerical recipes' random number generators.
00203 */
00204 
00205 #define NTAB 32                 /*   needs for ran1 & uniform_sample()   */
00206 #define EPS 1.2e-7              /*   needs for ran1 & uniform_sample()   */
00207 #define RNMX (1.0-EPS)          /*   needs for ran1 & uniform_sample()   */
00208 #define IM1 2147483563          /*   needs for uniform_sample()          */
00209 #define IM2 2147483399          /*   needs for uniform_sample()          */
00210 #define AM1 (1.0/IM1)           /*   needs for uniform_sample()          */
00211 #define IMM1 (IM1-1)            /*   needs for uniform_sample()          */
00212 #define IA1 40014               /*   needs for uniform_sample()          */
00213 #define IA2 40692               /*   needs for uniform_sample()          */
00214 #define IQ1 53668               /*   needs for uniform_sample()          */
00215 #define IQ2 52774               /*   needs for uniform_sample()          */
00216 #define IR1 12211               /*   needs for uniform_sample()          */
00217 #define IR2 3791                /*   needs for uniform_sample()          */
00218 #define NDIV1 (1+IMM1/NTAB)     /*   needs for uniform_sample()          */
00219 
00220 real uniform_sample()  
00221 {
00222     int j;
00223     int32_t k;
00224     static int32_t idum2=123456789;
00225     static int32_t iy=0;
00226     static int32_t iv[NTAB];
00227     real temp;
00228 
00229     if (the_seed <= 0) {
00230         if (-the_seed < 1) the_seed=1;
00231         else the_seed = -the_seed;
00232         idum2=the_seed;
00233         for (j=NTAB+7;j>=0;j--) {
00234             k=the_seed/IQ1;
00235             the_seed=IA1*(the_seed-k*IQ1)-k*IR1;
00236             if (the_seed < 0) the_seed += IM1;
00237             if (j < NTAB) iv[j] = the_seed;
00238         }
00239         iy=iv[0];
00240     }
00241     k=the_seed/IQ1;
00242     the_seed=IA1*(the_seed-k*IQ1)-k*IR1;
00243     if (the_seed < 0) the_seed += IM1;
00244     k=idum2/IQ2;
00245     idum2=IA2*(idum2-k*IQ2)-k*IR2;
00246     if (idum2 < 0) idum2 += IM2;
00247     j=int(iy/NDIV1);
00248     iy=iv[j]-idum2;
00249     iv[j] = the_seed;
00250     if (iy < 1) iy += IMM1;
00251     if ((temp=AM1*iy) > RNMX) return RNMX;
00252     else return temp;
00253 }
00254 
00255 /*  
00256     bounded_uniform(): return an uniform random generator in [a,b].
00257 */
00258 
00259 real  bounded_uniform(real a,real b)
00260 {
00261     real res = uniform_sample()*(b-a) + a;
00262     if (res >= b) return b*RNMX;
00263     else return res;
00264 }
00265 
00266 #undef NDIV
00267 #undef EPS
00268 #undef RNMX
00269 #undef IM1
00270 #undef IM2
00271 #undef AM1
00272 #undef IMM1
00273 #undef IA1
00274 #undef IA2
00275 #undef IQ1
00276 #undef IQ2
00277 #undef IR1
00278 #undef IR2
00279 #undef NDIV1
00280 
00281 /*  
00282     TRANSFORMATION METHOD:
00283     ---------------------
00284 */
00285 
00286 /*  
00287     expdev(): exponential deviate random number from the 'numerical recipes'.
00288 */
00289 
00290 real  expdev()
00291 {
00292     real dum;
00293 
00294     do
00295         dum=uniform_sample();
00296     while (fast_exact_is_equal(dum, 0.0));
00297     return -pl_log(dum);
00298 }
00299 
00300 real gaussian_01() 
00301 {
00302     real fac,rsq,v1,v2;
00303 
00304     if(the_seed < 0) iset=0;
00305     if (iset == 0) {
00306         do {
00307             v1=2.0*uniform_sample()-1.0;
00308             v2=2.0*uniform_sample()-1.0;
00309             rsq=v1*v1+v2*v2;
00310         } while (rsq >= 1.0 || fast_exact_is_equal(rsq, 0.0));
00311         fac=sqrt(-2.0*pl_log(rsq)/rsq);
00312         gset=v1*fac;
00313         iset=1;
00314         return v2*fac;
00315     } else {
00316         iset=0;
00317         return gset;
00318     }
00319 }
00320 
00321 /*  
00322     gaussian_mu_sigma(): returns a gaussian distributed random number
00323     with mean 'mu' and standard deviation 'sigma'.
00324 
00325     Rem: - i.e. N(mu,sigma).
00326 */
00327 
00328 real  gaussian_mu_sigma(real mu, real sigma)
00329 {
00330     return gaussian_01() * sigma + mu;
00331 }
00332 
00333 
00334 /*  
00335     gaussian_misture_mu_sigma(): returns a random number with mixture of gaussians,
00336     'w' is the mixture (must be positive summing to 1).
00337 
00338     Rem: - i.e. SUM w[i] * N(mu[i],sigma[i])
00339 */
00340 
00341 real  gaussian_mixture_mu_sigma(Vec& w, const Vec& mu, const Vec& sigma)
00342 {
00343     int    i;
00344     int    n = w.length();
00345     real *p_mu = mu.data();
00346     real *p_sigma = sigma.data();
00347     real *p_w = w.data();
00348     real  res = 0.0;
00349 
00350     for (i=0; i<n; i++, p_mu++, p_sigma++, p_w++)
00351         res += *p_w * gaussian_mu_sigma(*p_mu,*p_sigma);
00352 
00353     return res;
00354 }
00355 
00356 /*  
00357     REJECTION METHOD:
00358     ----------------
00359 */
00360 
00361 real  gamdev(int ia)
00362 {
00363     int j;
00364     real am,e,s,v1,v2,x,y;
00365 
00366     if (ia < 1) PLERROR("Error in routine gamdev");
00367     if (ia < 6) {
00368         x=1.0;
00369         for (j=1;j<=ia;j++) x *= uniform_sample();
00370         x = -pl_log(x);
00371     } else {
00372         do {
00373             do {
00374                 do {
00375                     v1=uniform_sample();
00376                     v2=2.0*uniform_sample()-1.0;
00377                 } while (v1*v1+v2*v2 > 1.0);
00378                 y=v2/v1;
00379                 am=ia-1;
00380                 s=sqrt(2.0*am+1.0);
00381                 x=s*y+am;
00382             } while (x <= 0.0);
00383             e=(1.0+y*y)*exp(am*pl_log(x/am)-s*y);
00384         } while (uniform_sample() > e);
00385     }
00386     return x;
00387 }
00388 
00389 real  poidev(real xm)
00390 {
00391     static real sq,alxm,g,oldm=(-1.0);
00392     real em,t,y;
00393 
00394     if (xm < 12.0) {
00395         if (!fast_exact_is_equal(xm, oldm)) {
00396             oldm=xm;
00397             g=exp(-xm);
00398         }
00399         em = -1;
00400         t=1.0;
00401         do {
00402             ++em;
00403             t *= uniform_sample();
00404         } while (t > g);
00405     } else {
00406         if (!fast_exact_is_equal(xm, oldm)) {
00407             oldm=xm;
00408             sq=sqrt(2.0*xm);
00409             alxm=pl_log(xm);
00410             g=xm*alxm-log_gamma(xm+1.0);
00411         }
00412         do {
00413             do {
00414                 y=tan(Pi*uniform_sample());
00415                 em=sq*y+xm;
00416             } while (em < 0.0);
00417             em=floor(em);
00418             t=0.9*(1.0+y*y)*exp(em*alxm-log_gamma(em+1.0)-g);
00419         } while (uniform_sample() > t);
00420     }
00421     return em;
00422 }
00423 
00424 /*  
00425     bnldev(): return a random deviate drawn from a binomial distribution of 'n' trials
00426     each of probability 'pp', from 'numerical recipes'.
00427 
00428     Rem: - The returned type is an real although a binomial random variable is an integer.
00429 */
00430 
00431 real  bnldev(real pp, int n)
00432 {
00433     int j;
00434     static int nold=(-1);
00435     real am,em,g,angle,p,bnl,sq,t,y;
00436     static real pold=(-1.0),pc,plog,pclog,en,oldg;
00437 
00438     p=(pp <= 0.5 ? pp : 1.0-pp);
00439     am=n*p;
00440     if (n < 25) {
00441         bnl=0.0;
00442         for (j=1;j<=n;j++)
00443             if (uniform_sample() < p) ++bnl;
00444     } else if (am < 1.0) {
00445         g=exp(-am);
00446         t=1.0;
00447         for (j=0;j<=n;j++) {
00448             t *= uniform_sample();
00449             if (t < g) break;
00450         }
00451         bnl=(j <= n ? j : n);
00452     } else {
00453         if (n != nold) {
00454             en=n;
00455             oldg=log_gamma(en+1.0);
00456             nold=n;
00457         } if (!fast_exact_is_equal(p, pold)) {
00458             pc=1.0-p;
00459             plog=pl_log(p);
00460             pclog=pl_log(pc);
00461             pold=p;
00462         }
00463         sq=sqrt(2.0*am*pc);
00464         do {
00465             do {
00466                 angle=Pi*uniform_sample();
00467                 y=tan(angle);
00468                 em=sq*y+am;
00469             } while (em < 0.0 || em >= (en+1.0));
00470             em=floor(em);
00471             t=1.2*sq*(1.0+y*y)*exp(oldg-log_gamma(em+1.0)
00472                                    -log_gamma(en-em+1.0)+em*plog+(en-em)*pclog);
00473         } while (uniform_sample() > t);
00474         bnl=em;
00475     }
00476     if (!fast_exact_is_equal(p, pp)) bnl=n-bnl;
00477     return bnl;
00478 }
00479 
00480 /*  
00481     SOME KIND OF DISCRETE DISTRIBUTIONS:
00482     -----------------------------------
00483 */
00484 
00485 /*  
00486     multinomial_sample(): returns a random deviate from a discrete distribution
00487     given explicitely by 'distribution'.
00488 
00489     Rem: - So, the vector elements of 'distribution' are probabilities summing to 1.
00490 
00491     - The returned value is a index value of 'distribution' (i.e. in range
00492     [0 .. (distribution->lenght)-1] ).
00493 
00494     - The graphical representation of vectors 'distribution' is histogram.
00495 
00496     - This random deviate is computed by the transformation method.
00497 */
00498 
00499 int  multinomial_sample(const Vec& distribution)
00500 {
00501     real  u  = uniform_sample();
00502     real* pi = distribution.data();
00503     real  s  = *pi;
00504     int    n  = distribution.length();
00505     int    i  = 0;
00506     while ((i<n) && (s<u)) {
00507         i++;
00508         pi++;
00509         s += *pi;
00510     }
00511     if (i==n)
00512         i = n - 1; /*   improbable but...   */
00513     return i;
00514 }
00515 
00516 int uniform_multinomial_sample(int N)
00517 {
00518     // N.B. uniform_sample() cannot return 1.0
00519     return int(N*uniform_sample());
00520 }
00521 
00522 void fill_random_uniform(const Vec& dest, real minval, real maxval)
00523 {
00524     Vec::iterator it = dest.begin();
00525     Vec::iterator itend = dest.end();  
00526     double scale = maxval-minval;
00527     for(; it!=itend; ++it)
00528         *it = real(uniform_sample()*scale+minval);
00529 }
00530 
00531 void fill_random_discrete(const Vec& dest, const Vec& set)
00532 {
00533     Vec::iterator it = dest.begin();
00534     Vec::iterator itend = dest.end();  
00535     int n=set.length();
00536     for(; it!=itend; ++it)
00537         *it = set[uniform_multinomial_sample(n)];
00538 }
00539 
00540 void fill_random_normal(const Vec& dest, real mean, real stdev)
00541 {
00542     Vec::iterator it = dest.begin();
00543     Vec::iterator itend = dest.end();  
00544     for(; it!=itend; ++it)
00545         *it = real(gaussian_mu_sigma(mean,stdev));
00546 }
00547 
00548 void fill_random_normal(const Vec& dest, const Vec& mean, const Vec& stdev)
00549 {
00550 #ifdef BOUNDCHECK
00551     if(mean.length()!=dest.length() || stdev.length()!=dest.length())
00552         PLERROR("In fill_random_normal: dest, mean and stdev must have the same length");
00553 #endif
00554     Vec::iterator it_mean = mean.begin();
00555     Vec::iterator it_stdev = stdev.begin();
00556     Vec::iterator it = dest.begin();
00557     Vec::iterator itend = dest.end();  
00558     for(; it!=itend; ++it, ++it_mean, ++it_stdev)
00559         *it = real(gaussian_mu_sigma(*it_mean,*it_stdev));
00560 }
00561 
00562 
00563 void fill_random_uniform(const Mat& dest, real minval, real maxval)
00564 { 
00565     double scale = maxval-minval;
00566     Mat::iterator it = dest.begin();
00567     Mat::iterator itend = dest.end();
00568     for(; it!=itend; ++it)
00569         *it = real(uniform_sample()*scale+minval); 
00570 }
00571 
00572 void fill_random_normal(const Mat& dest, real mean, real sdev)
00573 { 
00574     Mat::iterator it = dest.begin();
00575     Mat::iterator itend = dest.end();
00576     for(; it!=itend; ++it)
00577         *it = real(gaussian_mu_sigma(mean,sdev));
00578 }
00579 
00580 /*                                                      incbet.c
00581  *
00582  *      Incomplete beta integral
00583  *
00584  *
00585  * SYNOPSIS:
00586  *
00587  * double a, b, x, y, incbet();
00588  *
00589  * y = incbet( a, b, x );
00590  *
00591  *
00592  * DESCRIPTION:
00593  *
00594  * Returns incomplete beta integral of the arguments, evaluated
00595  * from zero to x.  The function is defined as
00596  *
00597  *                  x
00598  *     -            -
00599  *    | (a+b)      | |  a-1     b-1
00600  *  -----------    |   t   (1-t)   dt.
00601  *   -     -     | |
00602  *  | (a) | (b)   -
00603  *                 0
00604  *
00605  * The domain of definition is 0 <= x <= 1.  In this
00606  * implementation a and b are restricted to positive values.
00607  * The integral from x to 1 may be obtained by the symmetry
00608  * relation
00609  *
00610  *    1 - incbet( a, b, x )  =  incbet( b, a, 1-x ).
00611  *
00612  * The integral is evaluated by a continued fraction expansion
00613  * or, when b*x is small, by a power series.
00614  *
00615  * ACCURACY:
00616  *
00617  * Tested at uniformly distributed random points (a,b,x) with a and b
00618  * in "domain" and x between 0 and 1.
00619  *                                        Relative error
00620  * arithmetic   domain     # trials      peak         rms
00621  *    IEEE      0,5         10000       6.9e-15     4.5e-16
00622  *    IEEE      0,85       250000       2.2e-13     1.7e-14
00623  *    IEEE      0,1000      30000       5.3e-12     6.3e-13
00624  *    IEEE      0,10000    250000       9.3e-11     7.1e-12
00625  *    IEEE      0,100000    10000       8.7e-10     4.8e-11
00626  * Outputs smaller than the IEEE gradual underflow threshold
00627  * were excluded from these statistics.
00628  *
00629  * ERROR MESSAGES:
00630  *   message         condition      value returned
00631  * incbet domain      x<0, x>1          0.0
00632  * incbet underflow                     0.0
00633  */
00634 
00635 
00636 /*
00637   Cephes Math Library, Release 2.8:  June, 2000
00638   Copyright 1984, 1995, 2000 by Stephen L. Moshier
00639 */
00640 
00641 //#include "mconf.h"
00642 
00643 #define MAXGAM 171.624376956302725
00644 
00645 /*
00646   extern double MACHEP, MINLOG, MAXLOG;
00647   #ifdef ANSIPROT
00648   extern double gamma ( double );
00649   extern double lgam ( double );
00650   extern double exp ( double );
00651   extern double log ( double );
00652   extern double pow ( double, double );
00653   extern double fabs ( double );
00654   static double incbcf(double, double, double);
00655   static double incbd(double, double, double);
00656   static double pseries(double, double, double);
00657   #else
00658   double gamma(), lgam(), exp(), log(), pow(), fabs();
00659   static double incbcf(), incbd(), pseries();
00660   #endif
00661 
00662 */
00663 double MAXLOG =  7.09782712893383996732E2;     /* log(MAXNUM) */
00664 double MINLOG = -7.451332191019412076235E2;     /* log(2**-1075) */
00665 double MACHEP =  1.11022302462515654042E-16;   /* 2**-53 */
00666 //double pseries( double a, double b, double x );
00667 // double incbcf( double a, double b, double x ); // Does not seem to be used anymore.
00668 // double incbd( double a, double b, double x );  // Does not seem to be used anymore.
00669 double big = 4.503599627370496e15;
00670 double biginv =  2.22044604925031308085e-16;
00671 
00672 
00673 /*
00674   double incbet(double aa, double bb, double xx )
00675   {
00676   double a, b, t, x, xc, w, y;
00677   int flag;
00678 
00679   if( aa <= 0.0 || bb <= 0.0 )
00680   goto domerr;
00681 
00682   if( (xx <= 0.0) || ( xx >= 1.0) )
00683   {
00684   if( xx == 0.0 )
00685   return(0.0);
00686   if( xx == 1.0 )
00687   return( 1.0 );
00688   domerr:
00689   PLERROR("incbet: arguments out of expected domain");
00690   return( 0.0 );
00691   }
00692 
00693   flag = 0;
00694   if( (bb * xx) <= 1.0 && xx <= 0.95)
00695   {
00696   t = pseries(aa, bb, xx);
00697   goto done;
00698   }
00699 
00700   w = 1.0 - xx;
00701 
00702 // Reverse a and b if x is greater than the mean.
00703 if( xx > (aa/(aa+bb)) )
00704 {
00705 flag = 1;
00706 a = bb;
00707 b = aa;
00708 xc = xx;
00709 x = w;
00710 }
00711 else
00712 {
00713     a = aa;
00714     b = bb;
00715     xc = w;
00716     x = xx;
00717 }
00718 
00719 if( flag == 1 && (b * x) <= 1.0 && x <= 0.95)
00720 {
00721     t = pseries(a, b, x);
00722     goto done;
00723 }
00724 
00725 // Choose expansion for better convergence.
00726 y = x * (a+b-2.0) - (a-1.0);
00727 if( y < 0.0 )
00728     w = incbcf( a, b, x );
00729 else
00730 w = incbd( a, b, x ) / xc;
00731 
00732 // Multiply w by the factor
00733 //   a      b   _             _     _
00734 //  x  (1-x)   | (a+b) / ( a | (a) | (b) ) .
00735 
00736 y = a * log(x);
00737 t = b * log(xc);
00738 if( (a+b) < MAXGAM && fabs(y) < MAXLOG && fabs(t) < MAXLOG )
00739 {
00740     t = pow(xc,b);
00741     t *= pow(x,a);
00742     t /= a;
00743     t *= w;
00744     t *= gamma(a+b) / (gamma(a) * gamma(b));
00745     goto done;
00746 }
00747 // Resort to logarithms.
00748 y += t + log_gamma(a+b) - log_gamma(a) - log_gamma(b);
00749 y += log(w/a);
00750 if( y < MINLOG )
00751     t = 0.0;
00752 else
00753 t = exp(y);
00754 
00755 done:
00756 
00757 if( flag == 1 )
00758 {
00759     if( t <= MACHEP )
00760         t = 1.0 - MACHEP;
00761     else
00762         t = 1.0 - t;
00763 }
00764 return( t );
00765 }
00766 */
00767 
00768 /* Continued fraction expansion #1
00769  * for incomplete beta integral
00770  */
00771 
00772   double incbcf( double a, double b, double x )
00773 {
00774     double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
00775     double k1, k2, k3, k4, k5, k6, k7, k8;
00776     double r, t, ans, thresh;
00777     int n;
00778 
00779     k1 = a;
00780     k2 = a + b;
00781     k3 = a;
00782     k4 = a + 1.0;
00783     k5 = 1.0;
00784     k6 = b - 1.0;
00785     k7 = k4;
00786     k8 = a + 2.0;
00787 
00788     pkm2 = 0.0;
00789     qkm2 = 1.0;
00790     pkm1 = 1.0;
00791     qkm1 = 1.0;
00792     ans = 1.0;
00793     r = 1.0;
00794     n = 0;
00795     thresh = 3.0 * MACHEP;
00796     do
00797     {
00798         
00799         xk = -( x * k1 * k2 )/( k3 * k4 );
00800         pk = pkm1 +  pkm2 * xk;
00801         qk = qkm1 +  qkm2 * xk;
00802         pkm2 = pkm1;
00803         pkm1 = pk;
00804         qkm2 = qkm1;
00805         qkm1 = qk;
00806 
00807         xk = ( x * k5 * k6 )/( k7 * k8 );
00808         pk = pkm1 +  pkm2 * xk;
00809         qk = qkm1 +  qkm2 * xk;
00810         pkm2 = pkm1;
00811         pkm1 = pk;
00812         qkm2 = qkm1;
00813         qkm1 = qk;
00814 
00815         if( !fast_exact_is_equal(qk, 0) )
00816             r = pk/qk;
00817         if( !fast_exact_is_equal(r, 0) )
00818         {
00819             t = fabs( (ans - r)/r );
00820             ans = r;
00821         }
00822         else
00823             t = 1.0;
00824 
00825         if( t < thresh )
00826             goto cdone;
00827 
00828         k1 += 1.0;
00829         k2 += 1.0;
00830         k3 += 2.0;
00831         k4 += 2.0;
00832         k5 += 1.0;
00833         k6 -= 1.0;
00834         k7 += 2.0;
00835         k8 += 2.0;
00836 
00837         if( (fabs(qk) + fabs(pk)) > big )
00838         {
00839             pkm2 *= biginv;
00840             pkm1 *= biginv;
00841             qkm2 *= biginv;
00842             qkm1 *= biginv;
00843         }
00844         if( (fabs(qk) < biginv) || (fabs(pk) < biginv) )
00845         {
00846             pkm2 *= big;
00847             pkm1 *= big;
00848             qkm2 *= big;
00849             qkm1 *= big;
00850         }
00851     }
00852     while( ++n < 300 );
00853 
00854  cdone:
00855     return(ans);
00856 }
00857 
00858 /* Continued fraction expansion #2
00859  * for incomplete beta integral
00860  */
00861 
00862 /* Does not seem to be used anymore.
00863    double incbd( double a, double b, double x )
00864    {
00865    double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
00866    double k1, k2, k3, k4, k5, k6, k7, k8;
00867    double r, t, ans, z, thresh;
00868    int n;
00869 
00870    k1 = a;
00871    k2 = b - 1.0;
00872    k3 = a;
00873    k4 = a + 1.0;
00874    k5 = 1.0;
00875    k6 = a + b;
00876    k7 = a + 1.0;;
00877    k8 = a + 2.0;
00878 
00879    pkm2 = 0.0;
00880    qkm2 = 1.0;
00881    pkm1 = 1.0;
00882    qkm1 = 1.0;
00883    z = x / (1.0-x);
00884    ans = 1.0;
00885    r = 1.0;
00886    n = 0;
00887    thresh = 3.0 * MACHEP;
00888    do
00889    {
00890         
00891    xk = -( z * k1 * k2 )/( k3 * k4 );
00892    pk = pkm1 +  pkm2 * xk;
00893    qk = qkm1 +  qkm2 * xk;
00894    pkm2 = pkm1;
00895    pkm1 = pk;
00896    qkm2 = qkm1;
00897    qkm1 = qk;
00898 
00899    xk = ( z * k5 * k6 )/( k7 * k8 );
00900    pk = pkm1 +  pkm2 * xk;
00901    qk = qkm1 +  qkm2 * xk;
00902    pkm2 = pkm1;
00903    pkm1 = pk;
00904    qkm2 = qkm1;
00905    qkm1 = qk;
00906 
00907    if( qk != 0 )
00908    r = pk/qk;
00909    if( r != 0 )
00910    {
00911    t = fabs( (ans - r)/r );
00912    ans = r;
00913    }
00914    else
00915    t = 1.0;
00916 
00917    if( t < thresh )
00918    goto cdone;
00919 
00920    k1 += 1.0;
00921    k2 -= 1.0;
00922    k3 += 2.0;
00923    k4 += 2.0;
00924    k5 += 1.0;
00925    k6 += 1.0;
00926    k7 += 2.0;
00927    k8 += 2.0;
00928 
00929    if( (fabs(qk) + fabs(pk)) > big )
00930    {
00931    pkm2 *= biginv;
00932    pkm1 *= biginv;
00933    qkm2 *= biginv;
00934    qkm1 *= biginv;
00935    }
00936    if( (fabs(qk) < biginv) || (fabs(pk) < biginv) )
00937    {
00938    pkm2 *= big;
00939    pkm1 *= big;
00940    qkm2 *= big;
00941    qkm1 *= big;
00942    }
00943    }
00944    while( ++n < 300 );
00945    cdone:
00946    return(ans);
00947    }
00948 */
00949 
00950 /* Power series for incomplete beta integral.
00951    Use when b*x is small and x not too close to 1.  */
00952 
00953 /*
00954   double pseries( double a, double b, double x )
00955   {
00956   double s, t, u, v, n, t1, z, ai;
00957 
00958   ai = 1.0 / a;
00959   u = (1.0 - b) * x;
00960   v = u / (a + 1.0);
00961   t1 = v;
00962   t = u;
00963   n = 2.0;
00964   s = 0.0;
00965   z = MACHEP * ai;
00966   while( fabs(v) > z )
00967   {
00968   u = (n - b) * x / n;
00969   t *= u;
00970   v = t / (a + n);
00971   s += v; 
00972   n += 1.0;
00973   }
00974   s += t1;
00975   s += ai;
00976 
00977   u = a * log(x);
00978   if( (a+b) < MAXGAM && fabs(u) < MAXLOG )
00979   {
00980   t = gamma(a+b)/(gamma(a)*gamma(b));
00981   s = s * t * pow(x,a);
00982   }
00983   else
00984   {
00985   t = log_gamma(a+b) - log_gamma(a) - log_gamma(b) + u + log(s);
00986   if( t < MINLOG )
00987   s = 0.0;
00988   else
00989   s = exp(t);
00990   }
00991   return(s);
00992   }
00993 */
00994 
00995 void random_subset_indices(const TVec<int>& dest, int n)
00996 {
00997     if (dest.length()>n)
00998         PLERROR("random_subset_indices: 1st argument should have length (%d) <= value of 2nd argument (%d)",
00999                 dest.length(),n);
01000     TVec<int> v(0, n-1, 1);
01001     shuffleElements(v);
01002     dest << v.subVec(0,dest.length());
01003 }
01004 
01005 } // end of namespace PLearn
01006 
01007 
01008 /*
01009   Local Variables:
01010   mode:c++
01011   c-basic-offset:4
01012   c-file-style:"stroustrup"
01013   c-file-offsets:((innamespace . 0)(inline-open . 0))
01014   indent-tabs-mode:nil
01015   fill-column:79
01016   End:
01017 */
01018 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines