/* *
 * Title:        

 * Description:  

 * Copyright:    Copyright (c) 

 * Company:      

 * @author
 * @version 1.0
 */

package umontreal.iro.lecuyer.probdist;
import umontreal.iro.lecuyer.util.Num;


/**
 * Extends the class {@link ContinuousDistribution} for
 * the Johnson SB distribution
 * with shape parameters γ and 
 * δ > 0, location parameter ξ, 
 * and scale parameter λ > 0.
 * Denoting 
 * y = (x - ξ)/λ, the density is
 * 

 * 
 * 

 * f (x) = δ/(λy(1-y)(2π)1/2exp(- (1/2)[γ + δln(y/(1 - y))]2)) for ξ < x < ξ + λ,
 * 

 * and 0 elsewhere.  The distribution function is
 * 

 * 
 * 

 * F(x) = Φ[γ + δln(y/(1 - y))], for ξ < x < ξ + λ,
 * 

 * where Φ is the standard normal distribution function.
 * The inverse distribution function is
 * 

 * 
 * 

 * F-1(u) = ξ + λ(1/(1 + e-v(u)))        for 0 <= u <= 1,
 * 

 * where
 * 

 * 
 * 

 * v(u) = [Φ-1(u) - γ]/δ.
 * 

 * 
 * 

 * This class relies on the methods  {@link NormalDist#cdf01 NormalDist.cdf01} and
 *   {@link NormalDist#inverseF01 NormalDist.inverseF01}
 * of {@link NormalDist} to approximate Φ and Φ-1.
 * 
 */
public class JohnsonSBDist extends ContinuousDistribution {
   private double gamma;
   private double delta;
   private double xi;
   private double lambda;


   /**
    * Constructs a JohnsonSBDist object
    *    with shape parameters γ and δ,
    *    location parameter ξ and scale parameter λ.
    * 
    */
   public JohnsonSBDist (double gamma, double delta,
                         double xi, double lambda) {
      setParams (gamma, delta, xi, lambda);
   }


   public double density (double x) {
      return density (gamma, delta, xi, lambda, x);
   }
       
   public double cdf (double x) {
      return cdf (gamma, delta, xi, lambda, x);
   }
     
   public double barF (double x) {
      return barF (gamma, delta, xi, lambda, x);
   }

   public double inverseF (double u){
      return inverseF (gamma, delta, xi, lambda, u);
   }


   /**
    * Computes the density function.
    * 
    */
   public static double density (double gamma, double delta,
                                 double xi, double lambda, double x) {
      if (lambda <= 0)
         throw new IllegalArgumentException ("lambda <= 0");
      if (delta <= 0)
         throw new IllegalArgumentException ("delta <= 0");
      if (x <= xi || x >= (xi+lambda))
         return 0.0;
      double y = (x - xi)/lambda;
      double t1 = gamma + delta*Math.log (y/(1.0 - y));
      return delta/(lambda*y*(1.0 - y)*Math.sqrt (2.0*Math.PI))*
           Math.exp (-t1*t1/2.0);
   }
      

   /**
    * Computes the distribution function.
    * 
    */
   public static double cdf (double gamma, double delta,
                             double xi, double lambda, double x) {
      if (lambda <= 0)
         throw new IllegalArgumentException ("lambda <= 0");
      if (delta <= 0)
         throw new IllegalArgumentException ("delta <= 0");
      if (x <= xi)
         return 0.0;
      if (x >= xi+lambda)
         return 1.0;
      double y = (x - xi)/lambda;
      return NormalDist.cdf01 (gamma + delta*Math.log (y/(1.0 - y)));
   }


   /**
    * Computes the complementary distribution.
    * 
    */
   public static double barF (double gamma, double delta,
                              double xi, double lambda, double x) {
      if (lambda <= 0)
         throw new IllegalArgumentException ("lambda <= 0");
      if (delta <= 0)
         throw new IllegalArgumentException ("delta <= 0");
      if (x <= xi)
         return 1.0;
      if (x >= xi+lambda)
         return 0.0;
      double y = (x - xi)/lambda;
      return NormalDist.barF01 (gamma + delta*Math.log (y/(1.0 - y)));
   }


   /**
    * Computes the inverse of the distribution.
    * 
    */
   public static double inverseF (double gamma, double delta,
                                  double xi, double lambda, double u) {
      if (lambda <= 0)
         throw new IllegalArgumentException ("lambda <= 0");
      if (delta <= 0)
         throw new IllegalArgumentException ("delta <= 0");
      if (u > 1.0 || u < 0.0)
          throw new IllegalArgumentException ("u not in [0,1]");

      double v, z;

      if (u >= 1.0)    // if u == 1, in fact
          return xi+lambda;
      if (u <= 0.0)    // if u == 0, in fact
          return xi;

      z = NormalDist.inverseF01 (u);
      v = (z - gamma)/delta;

      if ((z >= XBIG) || (v >= Num.DBL_MAX_EXP*Num.LN2))
            return xi + lambda;

        if ((z <= -XBIG) || (v <= -Num.DBL_MAX_EXP*Num.LN2))
            return xi;

        v = Math.exp (v);
        return (xi + (xi+lambda)*v)/(1.0 + v);
   }


   /**
    * Returns the value of γ for this object.
    * 
    */
   public double getGamma() {
      return gamma;
   }
      

   /**
    * Returns the value of δ for this object.
    * 
    */
   public double getDelta() {
      return delta;
   }
      

   /**
    * Returns the value of ξ for this object.
    * 
    */
   public double getXi() {
      return xi;
   }
      

   /**
    * Returns the value of λ for this object.
    * 
    */
   public double getLambda() {
      return lambda;
   }

      

   /**
    * Sets the value of the parameters γ, δ, ξ and
    *   λ for this object.
    * 
    */
   public void setParams(double gamma, double delta, double xi,
                         double lambda) {
      if (lambda <= 0)
         throw new IllegalArgumentException ("lambda <= 0");
      if (delta <= 0)
         throw new IllegalArgumentException ("delta <= 0");
      this.gamma = gamma;
      this.delta = delta;
      this.xi = xi;
      this.lambda = lambda;
      supportA = xi;
      supportB = xi + lambda;
   }


   /**
    * Return a table containing the parameters of the current distribution.
    *    This table is put in regular order: [γ, δ, ξ, λ].
    * 
    * 
    */
   public double[] getParams () {
      double[] retour = {gamma, delta, xi, lambda};
      return retour;
   }


   public String toString () {
      return getClass().getSimpleName() + " : gamma = " + gamma + ", delta = " + delta + ", xi = " + xi + ", lambda = " + lambda;
   }

}