umontreal.iro.lecuyer.randvar.KernelDensityVarCorrectGen Maven / Gradle / Ivy

Go to download

Show more of this group Show more artifacts with this name
Show all versions of ssj Show documentation

SSJ is a Java library for stochastic simulation, developed under the direction of Pierre L'Ecuyer, in the Département d'Informatique et de Recherche Opérationnelle (DIRO), at the Université de Montréal. It provides facilities for generating uniform and nonuniform random variates, computing different measures related to probability distributions, performing goodness-of-fit tests, applying quasi-Monte Carlo methods, collecting (elementary) statistics, and programming discrete-event simulations with both events and processes.

The newest version!



/*
 * Class:        KernelDensityVarCorrectGen
 * Description:  random variate generators for distributions obtained via
                 kernel density estimation methods
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Université de Montréal
 * Organization: DIRO, Université de Montréal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   GPL licence site.
 */

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 s_n² 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 
 * σ_e² = 1 + (hσ_k/s_n)².
 * The constant 
 * σ_k² 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 
    * σ_k² = 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;
   }
}