package umontreal.iro.lecuyer.randvar;
   import umontreal.iro.lecuyer.functions.MathFunction;
   import umontreal.iro.lecuyer.rng.RandomStream;
   import umontreal.iro.lecuyer.probdist.ContinuousDistribution;
   import umontreal.iro.lecuyer.probdist.InverseDistFromDensity;



/**
 * Implements a method for generating random variates by numerical inversion of
 * an arbitrary continuous distribution when only the probability density
 * is known.  The cumulative probabilities (cdf) are pre-computed by
 *  numerical quadrature  of the
 * density using Gauss-Lobatto integration over suitably small intervals to
 * satisfy the required precision, and these values are kept in tables. Then the
 *  algorithm uses polynomial interpolation  over the tabulated values to get
 *   the inverse cdf. The user can select the
 *   desired precision and the degree of the interpolating polynomials.
 * 
 * 

 * The algorithm may fail for some distributions for which the density
 *  becomes infinite at a point (for ex. the Gamma and the Beta distributions
 *  with 
 * α < 1)  if one requires too high a precision
 * (a  too small eps, for ex. 
 * ε∼10-15).
 * However, it should work also for continuous densities with finite discontinuities.
 * 
 * 

 * While the setup time is relatively slow, the generation of random variables
 * is extremely fast and practically independent of the required precision
 * and of the specific distribution. The following table shows the time needed
 * (in seconds) to generate 108 random numbers using inversion
 * from a given class, then the numerical inversion with Gauss-Lobatto integration
 * implemented here, and finally the speed ratios between the two methods. 
 * The speed ratio is  the speed of the latter over the former.
 * Thus for the beta distribution with parameters (5, 500), generating random
 * variables with the Gauss-Lobatto integration implemented in this class is 
 * more than 1700
 * times faster than using inversion from the BetaDist class.
 * These tests were made on a machine with processor AMD Athlon 4000, running
 * Red Hat Linux, with clock speed at 2403 MHz.
 * 
 * 

 * 

 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
Distribution
Inversion
Gauss-Lobatto
speed ratio
NormalDist(10.5, 5)
9.19
8.89
1.03
ExponentialDist(5)
17.72
8.82
2.0
CauchyDist(10.5, 5)
18.30
8.81
2.1
BetaSymmetricalDist(10.5)
242.80
8.85
27.4
GammaDist(55)
899.50
8.89
101
ChiSquareNoncentralDist(10.5, 5)
5326.90
8.85
602
BetaDist(5, 500)
15469.10
8.86
1746
 * 
 * 
 * 

 * The following table gives the time (in sec.) needed to
 * create an object (setup time) and to generate one random variable for 
 * this class compared to the same for the inversion method specific to each class,
 * and the ratios of the times (init + one random variable) of the two methods.
 * For inversion, we initialized 108 times; for this class, we 
 * initialized 104 times.
 * 
 * 

 * 

 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
Distribution
Inversion
Gauss-Lobatto
time ratio
 
108 init
104 init
for 1 init
NormalDist(10.5, 5)
5.30
38.29
26426
ExponentialDist(5)
3.98
27.05
12466
CauchyDist(10.5, 5)
5.05
58.39
25007
BetaSymmetricalDist(10.5)
90.66
68.33
2049
GammaDist(55)
13.15
58.34
639
ChiSquareNoncentralDist(10.5, 5)
190.48
248.98
451
BetaDist(5, 500)
63.60
116.57
75
 * 
 * 
 * 

 * If only a few random variables are needed, then using this class
 * is not efficient because of the slow set-up. But if one wants to generate
 * large samples from the same distribution with fixed parameters, then
 * this class will be very efficient. The following table gives the number of
 * random variables generated beyond which, using this class will be
 * worthwhile. 
 * 
 * 

 * 

 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
Distribution
number of generated variables
NormalDist(10.5, 5)
41665
ExponentialDist(5)
15266
CauchyDist(10.5, 5)
31907
BetaSymmetricalDist(10.5)
2814
GammaDist(55)
649
ChiSquareNoncentralDist(10.5, 5)
467
BetaDist(5, 500)
75
 * 
 * 
 * 

 * Thus, for example, if one needs to generate less than 15266 exponential
 * random variables, then using the InverseFromDensityGen class
 * is not wortwhile: it will be faster to use inversion from the
 * ExponentialGen class.
 * 
 */
public class InverseFromDensityGen extends RandomVariateGen  {


   /**
    * Creates a new generator for the continuous distribution
    * dis, using stream s. dis must have a well-defined
    * density method; its other methods are unused. For a non-standard distribution
    * dis, the user may wish to set the left and the right boundaries
    * between which the density is non-zero by calling methods
    * {@link umontreal.iro.lecuyer.probdist.ContinuousDistribution#setXinf setXinf}
    * and
    * {@link umontreal.iro.lecuyer.probdist.ContinuousDistribution#setXsup setXsup}
    * of dis, for better efficiency.
    * Argument xc can be the mean,
    * the mode or any other x for which the density is relatively large.
    * The u-resolution eps is the desired absolute error in the CDF,
    * and order is the degree of the
    * Newton interpolating polynomial over each interval.
    * An order of 3 or 5, and an eps of 10-6 to 10-12
    * are usually good choices.
    *  Restrictions: 
    * 3 <= order <= 12.
    * 
    */
   public InverseFromDensityGen (RandomStream s, ContinuousDistribution dis,
                                 double xc, double eps, int order)  {
      super (s, null);
      dist = new InverseDistFromDensity (dis, xc, eps, order);
   }


   /**
    * Creates a new generator from the continuous probability density
    * dens. The left and the right boundaries of the density are
    * xleft and xright (the density is 0 outside the
    * interval [xleft, xright]).
    * See the description of the other constructor.
    * 
    */
   public InverseFromDensityGen (RandomStream s, MathFunction dens,
                                 double xc, double eps, int order,
                                 double xleft, double xright)  {
      super (s, null);
      dist = new InverseDistFromDensity (dens, xc, eps, order, xleft, xright);
   } 


   /**
    * Generates a new random variate.
    * 
    */
   public double nextDouble()  {
      return dist.inverseF (stream.nextDouble());
   }


   /**
    * Returns the xc given in the constructor.
    * 
    */
   public double getXc() {
      return ((InverseDistFromDensity)dist).getXc();
   }


   /**
    * Returns the u-resolution eps associated with this object.
    * 
    */
   public double getEpsilon() {
      return ((InverseDistFromDensity)dist).getEpsilon();
   }



   /**
    * Returns the order associated with this object.
    * 
    * 
    */
   public int getOrder() {
      return ((InverseDistFromDensity)dist).getOrder();
   }


}