PLearn 0.1
|
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(<ime); 00188 today = localtime(<ime); 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 :