package umontreal.iro.lecuyer.randvar;
import umontreal.iro.lecuyer.probdist.*;
import umontreal.iro.lecuyer.rng.RandomStream;

/**
 * This class implements random variate generators for distributions
 * obtained via kernel density estimation methods from a set of
 * n individual observations 
 * x1,..., xn. 
 * The basic idea is to center a copy of the same symmetric density 
 * at each observation and take an equally weighted mixture of the n
 * copies as an estimator of the density from which the observations come.
 * The resulting kernel density has the general form
 * 

 * 
 * 

 * fn(x) = (1/nh)∑i=1nk((x - xi)/h),
 * 

 * where k is a fixed pre-selected density called the kernel 
 * and h is a positive constant called the bandwidth or
 * smoothing factor.
 * A difficult practical issue is the selection of k and h.
 * Several approaches have been proposed for that; see, e.g.,.
 * 
 * 

 * The constructor of a generator from a kernel density requires a
 * random stream s, the n observations in the form of an empirical
 * distribution, a random variate generator for the kernel density k,
 * and the value of the bandwidth h.
 * The random variates are then generated as follows:
 * select an observation xI at random, by inversion, using stream s,
 * then generate random variate Y with the generator provided for
 * the density k, and return xI + hY.
 * 
 * 

 * A simple formula for the bandwidth, suggested in,
 * is 
 * h = αkh0, where
 * 

 * 
 * 

 * h0 = 1.36374min(sn, q/1.34)n-1/5,
 * 

 * sn and q are the empirical standard deviation and the 
 * interquartile range of the n observations,
 * and αk is a constant that depends on the type of kernel k.
 * It is defined by
 * 

 * 
 * 

 * αk = (σk-4∫-∞∞k(x)dx)1/5
 * 

 * where σk is the standard deviation of the density k.
 * The static method {@link #getBaseBandwidth getBaseBandwidth} permits one to compute h0 
 * for a given empirical distribution.
 * 
 * 

 * 

 * 

 * 

 * 
 * 
 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
Table:
 * Some suggested kernels
name
constructor
αk

 * σk2
efficiency
Epanechnikov
BetaSymmetricalDist(2, -1, 1)
1.7188
1/5
1.000
triangular
TriangularDist(-1, 1, 0)
1.8882
1/6
0.986
Gaussian
NormalDist()
0.7764
1
0.951
boxcar
UniformDist(-1, 1)
1.3510
1/3
0.930
logistic
LogisticDist()
0.4340
3.2899
0.888
Student-t(3)
StudentDist(3)
0.4802
3
0.674
 * 
 * 

 * 
 * 

 * Table  gives the precomputed values of σk and
 * αk for selected (popular) kernels.
 * The values are taken from.
 * The second column gives the name of a function (in this package)
 * that constructs the corresponding distribution.
 * The efficiency of a kernel is defined as the ratio of its
 * mean integrated square error over that of the Epanechnikov kernel,
 * which has optimal efficiency and corresponds to the beta distribution
 * with parameters (2, 2) over the interval (- 1, 1).
 * 
 */
public class KernelDensityGen extends RandomVariateGen {

   protected RandomVariateGen kernelGen;
   protected double bandwidth;
   protected boolean positive;   // If we want positive reflection.



   /**
    * Creates a new generator for a kernel density estimated 
    *     from the observations given by the empirical distribution dist,
    *     using stream s to select the observations,
    *     generator kGen to generate the added noise from the kernel
    *     density, and bandwidth h.
    * 
    */
   public KernelDensityGen (RandomStream s, EmpiricalDist dist,
                            RandomVariateGen kGen, double h) {
      super (s, dist);
      if (h < 0.0)
         throw new IllegalArgumentException ("h < 0");
      if (kGen == null)
         throw new IllegalArgumentException ("kGen == null");
      kernelGen = kGen;
      bandwidth = h;
   }


   /**
    * This constructor uses a gaussian kernel and the default 
    *     bandwidth 
    * h = αkh0 with the αk 
    *     suggested in Table  for the gaussian distribution.
    *     This kernel has an efficiency of 0.951.
    * 
    */
   public KernelDensityGen (RandomStream s, EmpiricalDist dist,
                            NormalGen kGen) {
      this (s, dist, kGen, 0.77639 * getBaseBandwidth (dist));
   }


   /**
    * Computes and returns the value of h0 in.
    * 
    */
   public static double getBaseBandwidth (EmpiricalDist dist) {
      double r = dist.getInterQuartileRange() / 1.34;
      double sigma = dist.getSampleStandardDeviation();
      if (sigma < r) r = sigma;
      return (1.36374 * r / Math.exp (0.2 * Math.log (dist.getN())));
   }


   /**
    * Sets the bandwidth to h.
    * 
    */
   public void setBandwidth (double h) {
      if (h < 0)
         throw new IllegalArgumentException ("h < 0");
      bandwidth = h;
   }


   /**
    * After this method is called with true,
    *   the generator will produce only positive values, by using 
    *   the reflection method: replace all negative values by their
    *   absolute values.
    *   That is, {@link #nextDouble nextDouble} will return | x| if x is the 
    *   generated variate.  The mecanism is disabled when the method is
    *   called with false.
    * 
    */
   public void setPositiveReflection (boolean reflect) {
      positive = reflect;
   }


   public double nextDouble() {
      double x = (dist.inverseF (stream.nextDouble())
                  + bandwidth * kernelGen.nextDouble());
      if (positive)
         return Math.abs (x);
      else
         return x;
   }
}