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


/**
 * .
 * 
 * This class implements a  (NIG) random variate generator by
 * using a normal generator (N) and an inverse gaussian generator (IG),
 *  as described in the following
 * 

 * 

 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
Y
∼
IG(δ/γ, δ2)

 *  
X | (Y = y)
∼
N(μ + βy, y).

 *  
 * 

 * 
 * The normal 
 * N(μ, σ2) has mean μ and variance σ2, while
 * the inverse gaussian has the parametrization described in{@link umontreal.iro.lecuyer.randvar.InverseGaussianGen InverseGaussianGen}.
 * If 
 * γ = (α^2 - β^2)1/2 with 
 * 0 <= | β| < α and
 * δ > 0, then 
 * X∼NIG(α, β, μ, δ). 
 * 
 */
public class NormalInverseGaussianIGGen extends NormalInverseGaussianGen  {
   private NormalGen genN;
   private InverseGaussianGen genIG;


   /**
    * Creates a  random variate generator with parameters
    *   α, β = beta, μ = mu and δ,
    *  using generators  ig and ng, as described
    * above. The parameters α and δ are included in generator ig.
    * 
    */
   public NormalInverseGaussianIGGen (InverseGaussianGen ig, NormalGen ng,
                                      double beta, double mu)  {
      super (null, null);
      setParams (ig, ng, beta, mu);
   }


   /**
    * Generates a new variate from the  distribution with 
    * parameters α, β = beta, μ = mu and δ, 
    * using generators ig and ng, as described in eq..
    * The parameters α and δ are included in generator ig.
    * 
    */
   public static double nextDouble (InverseGaussianGen ig, NormalGen ng,
                                    double beta, double mu)  {
      return mynig (ig, ng, beta, mu);
   }
 

   public double nextDouble() {
      return mynig (genIG, genN, beta, mu);
   }


// >>>>>>>>>>>>>>>>>>>>  P R I V A T E     M E T H O D S   <<<<<<<<<<<<<<<<<<<

   private static double mynig (InverseGaussianGen ig, NormalGen ng,
                                double beta, double mu) {

      double y = ig.nextDouble ();
      double x = mu + beta*y + Math.sqrt(y)*ng.nextDouble ();
      return x;
   }


   protected void setParams (InverseGaussianGen ig, NormalGen ng,
                             double beta, double mu) {
      delta = Math.sqrt(ig.getLambda());
      gamma = delta / ig.getMu();
      alpha = Math.sqrt(gamma*gamma + beta*beta);
      setParams (alpha, beta, mu, delta);
      this.genN = ng;
      this.genIG = ig;
   }
}