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

/**
 * This class is a variant of {@link KernelDensityGen}, but with
 * a rescaling of the empirical distribution so that the variance
 * of the density used to generate the random variates is equal
 * to the empirical variance,
 *  as suggested by Silverman.
 * 
 * 

 * Let bar(x)n and sn2 be the sample mean and sample variance 
 * of the observations.
 * The distance between each generated random variate and the 
 * sample mean bar(x)n is multiplied by the correcting factor
 * 
 * 1/σe, where 
 * σe2 = 1 + (hσk/sn)2.
 * The constant 
 * σk2 must be passed to the constructor.
 * Its value can be found in 
 * the Table in {@link KernelDensityGen} for some popular
 * kernels.
 * 
 */
public class KernelDensityVarCorrectGen extends KernelDensityGen {

   protected double sigmak2;   // Value of sigma_k^2.
   protected double mean;      // Sample mean of the observations.
   protected double invSigmae; // 1 / sigma_e.



   /**
    * 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, bandwidth h, and 
    * σk2 = sigmak2 used for
    *     the variance correction.
    * 
    */
   public KernelDensityVarCorrectGen (RandomStream s, EmpiricalDist dist,
                                      RandomVariateGen kGen, double h,
                                      double sigmak2) {
      super (s, dist, kGen, h);
      this.sigmak2 = sigmak2;
      mean = dist.getSampleMean();
      double var = dist.getSampleVariance();
      invSigmae = 1.0 / Math.sqrt (1.0 + h * h * sigmak2 / var);
   }


   /**
    * This constructor uses a gaussian kernel and the default 
    *     bandwidth suggested in Table  for the gaussian 
    *     distribution.
    * 
    */
   public KernelDensityVarCorrectGen (RandomStream s, EmpiricalDist dist,
                                      NormalGen kGen) {
      this (s, dist, kGen, 0.77639 * getBaseBandwidth (dist), 1.0);
   }


   public void setBandwidth (double h) {
      if (h < 0)
         throw new IllegalArgumentException ("h < 0");
      bandwidth = h;
      double var = ((EmpiricalDist) dist).getSampleVariance();
      invSigmae = 1.0 / Math.sqrt (1.0 + h * h * sigmak2 / var);
   }

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