1.1 --- a/lemon/random.h Mon Oct 01 18:55:58 2007 +0000
1.2 +++ b/lemon/random.h Mon Oct 01 18:57:21 2007 +0000
1.3 @@ -748,6 +748,44 @@
1.4 return -std::log(real<double>())/lambda;
1.5 }
1.6
1.7 + double gamma(int k)
1.8 + {
1.9 + double s = 0;
1.10 + for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());
1.11 + return s;
1.12 + }
1.13 +
1.14 + /// Gamma distribution with given shape and scale parameter
1.15 +
1.16 + /// This function generates a gamma distribution random number.
1.17 + ///
1.18 + ///\param k shape parameter (<tt>k>0</tt>)
1.19 + ///\param theta scale parameter
1.20 + ///
1.21 + double gamma(double k,double theta=1.0)
1.22 + {
1.23 + double xi,nu;
1.24 + const double delta = k-std::floor(k);
1.25 + const double v0=M_E/(M_E-delta);
1.26 + do {
1.27 + double V0=1.0-real<double>();
1.28 + double V1=1.0-real<double>();
1.29 + double V2=1.0-real<double>();
1.30 + if(V2<=v0)
1.31 + {
1.32 + xi=std::pow(V1,1.0/delta);
1.33 + nu=V0*std::pow(xi,delta-1.0);
1.34 + }
1.35 + else
1.36 + {
1.37 + xi=1.0-std::log(V1);
1.38 + nu=V0*std::exp(-xi);
1.39 + }
1.40 + } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));
1.41 + return theta*(xi-gamma(int(std::floor(k))));
1.42 + }
1.43 +
1.44 +
1.45 ///@}
1.46
1.47 ///\name Two dimensional distributions