JohnsonSBDist.java

/*
 * Class:        JohnsonSBDist
 * Description:  Johnson S_B distribution
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author
 * @since
 *
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
package umontreal.ssj.probdist;
 
import umontreal.ssj.util.Num;

public class JohnsonSBDist extends JohnsonSystem {
   // m_psi is used in computing the mean and the variance
   private double m_psi = -1.0e100;
 
   private static double getMeanPsi(double gamma, double delta, double xi, double lambda, double[] tpsi) {
      // Returns the theoretical mean of t = (x - xi)/lambda;
      // also compute psi and returns it in tpsi[0], since
      // it is used in computing the mean and the variance
 
      final int NMAX = 10000;
      final double EPS = 1.0e-15;
 
      double a1 = 1.0 / (2 * delta * delta);
      double a2 = (1.0 - 2 * delta * gamma) / (2 * delta * delta);
      double a3 = (gamma - 1. / delta) / delta;
      int n = 0;
      double tem = 1;
      double sum = 0;
      double v;
      while (Math.abs(tem) > EPS * Math.abs(sum) && n < NMAX) {
         ++n;
         v = Math.exp(-n * gamma / delta) + Math.exp(n * a3);
         tem = Math.exp(-n * n * a1) * v / (1 + Math.exp(-2 * n * a1));
         // tem = Math.exp(-n*n*a1) * Math.cosh(n*a2) / Math.cosh(n*a1);
         sum += tem;
      }
      if (n >= NMAX)
         System.err.println("JohnsonSBDist:  possible lack of accuracy on mean");
      double A = (0.5 + sum) / (delta);
 
      a1 = Math.PI * Math.PI * delta * delta;
      a2 = Math.PI * delta * gamma;
      int j;
      n = 0;
      tem = 1;
      sum = 0;
      while (Math.abs(tem) > EPS * Math.abs(sum) && n < NMAX) {
         ++n;
         j = 2 * n - 1;
         tem = Math.exp(-j * j * a1 / 2.0) * Math.sin(j * a2) / Math.sinh(j * a1);
         sum += tem;
      }
      if (n >= NMAX)
         System.err.println("JohnsonSBDist:  possible lack of accuracy on mean");
      double B = 2.0 * Math.PI * delta * sum;
 
      a1 = 2 * Math.PI * Math.PI * delta * delta;
      a2 = 2 * Math.PI * delta * gamma;
      n = 0;
      tem = 1;
      sum = 0;
      while (Math.abs(tem) > EPS * Math.abs(sum) && n < NMAX) {
         ++n;
         tem = Math.exp(-n * n * a1) * Math.cos(n * a2);
         sum += tem;
      }
      if (n >= NMAX)
         System.err.println("JohnsonSBDist:  possible lack of accuracy on mean");
      double C = 1 + 2.0 * sum;
 
      double D = Math.sqrt(2 * Math.PI) * Math.exp(gamma * gamma / 2.0);
      double tmean = (A - B) / (C * D);
      tpsi[0] = C * D;
      return tmean;
   }

   public JohnsonSBDist(double gamma, double delta, double xi, double lambda) {
      super(gamma, delta, xi, lambda);
      setLastParams(xi, lambda);
   }
 
   private void setLastParams(double xi, double lambda) {
      supportA = xi;
      supportB = 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);
   }
 
   public double getMean() {
      return JohnsonSBDist.getMean(gamma, delta, xi, lambda);
   }
 
   public double getVariance() {
      return JohnsonSBDist.getVariance(gamma, delta, xi, lambda);
   }
 
   public double getStandardDeviation() {
      return JohnsonSBDist.getStandardDeviation(gamma, delta, xi, lambda);
   }

   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 z = gamma + delta * Math.log(y / (1.0 - y));
      return delta / (lambda * y * (1.0 - y) * Math.sqrt(2.0 * Math.PI)) * Math.exp(-z * z / 2.0);
   }

   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;
      double z = gamma + delta * Math.log(y / (1.0 - y));
      return NormalDist.cdf01(z);
   }

   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;
      double z = gamma + delta * Math.log(y / (1.0 - y));
      return NormalDist.barF01(z);
   }

   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]");
 
      if (u >= 1.0) // if u == 1, in fact
         return xi + lambda;
      if (u <= 0.0) // if u == 0, in fact
         return xi;
 
      double z = NormalDist.inverseF01(u);
      double v = (z - gamma) / delta;
 
      if (v >= Num.DBL_MAX_EXP * Num.LN2)
         return xi + lambda;
      if (v <= Num.DBL_MIN_EXP * Num.LN2)
         return xi;
 
      v = Math.exp(v);
      return (xi + (xi + lambda) * v) / (1.0 + v);
   }

   public static double[] getMLE(double[] x, int n, double xi, double lambda) {
      if (n <= 0)
         throw new IllegalArgumentException("n <= 0");
      final double LN_EPS = Num.LN_DBL_MIN - Num.LN2;
      double[] ftab = new double[n];
      double sum = 0.0;
      double t;
 
      for (int i = 0; i < n; i++) {
         t = (x[i] - xi) / lambda;
         if (t <= 0.)
            ftab[i] = LN_EPS;
         else if (t >= 1 - Num.DBL_EPSILON)
            ftab[i] = Math.log(1. / Num.DBL_EPSILON);
         else
            ftab[i] = Math.log(t / (1. - t));
         sum += ftab[i];
      }
      double empiricalMean = sum / n;
 
      sum = 0.0;
      for (int i = 0; i < n; i++) {
         t = ftab[i] - empiricalMean;
         sum += t * t;
      }
      double sigmaf = Math.sqrt(sum / n);
 
      double[] param = new double[2];
      param[0] = -empiricalMean / sigmaf;
      param[1] = 1.0 / sigmaf;
 
      return param;
   }

   public static JohnsonSBDist getInstanceFromMLE(double[] x, int n, double xi, double lambda) {
      double parameters[] = getMLE(x, n, xi, lambda);
      return new JohnsonSBDist(parameters[0], parameters[1], xi, lambda);
   }

   public static double getMean(double gamma, double delta, double xi, double lambda) {
      if (lambda <= 0)
         throw new IllegalArgumentException("lambda <= 0");
      if (delta <= 0)
         throw new IllegalArgumentException("delta <= 0");
      double[] tpsi = new double[1];
      double mu = getMeanPsi(gamma, delta, xi, lambda, tpsi);
      return xi + lambda * mu;
   }

   public static double getVariance(double gamma, double delta, double xi, double lambda) {
      if (lambda <= 0.0)
         throw new IllegalArgumentException("lambda <= 0");
      if (delta <= 0.0)
         throw new IllegalArgumentException("delta <= 0");
 
      final int NMAX = 10000;
      final double EPS = 1.0e-15;
 
      double a1 = 1.0 / (2.0 * delta * delta);
      double a2 = (1.0 - 2.0 * delta * gamma) / (2.0 * delta * delta);
      double a3 = (gamma - 1. / delta) / delta;
      double v;
      int n = 0;
      double tem = 1;
      double sum = 0;
      while (Math.abs(tem) > EPS * Math.abs(sum) && n < NMAX) {
         ++n;
         v = Math.exp(-n * gamma / delta) - Math.exp(n * a3);
         tem = n * Math.exp(-n * n * a1) * v / (1 + Math.exp(-2 * n * a1));
         // tem = n*Math.exp(-n*n*a1) * Math.sinh(n*a2) / Math.cosh(n*a1);
         sum += tem;
      }
      if (n >= NMAX)
         System.err.println("JohnsonSBDist:  possible lack of accuracy on variance");
      double A = -sum / (delta * delta);
 
      a1 = Math.PI * Math.PI * delta * delta;
      a2 = Math.PI * delta * gamma;
      int j;
      n = 0;
      tem = 1;
      sum = 0;
      while (Math.abs(tem) > EPS * Math.abs(sum) && n < NMAX) {
         ++n;
         j = 2 * n - 1;
         tem = j * Math.exp(-j * j * a1 / 2.0) * Math.cos(j * a2) / Math.sinh(j * a1);
         sum += tem;
      }
      if (n >= NMAX)
         System.err.println("JohnsonSBDist:  possible lack of accuracy on variance");
      double B = 2.0 * a1 * sum;
 
      a1 = 2 * Math.PI * Math.PI * delta * delta;
      a2 = 2 * Math.PI * delta * gamma;
      n = 0;
      tem = 1;
      sum = 0;
      while (Math.abs(tem) > EPS * Math.abs(sum) && n < NMAX) {
         ++n;
         tem = n * Math.exp(-n * n * a1) * Math.sin(n * a2);
         sum += tem;
      }
      if (n >= NMAX)
         System.err.println("JohnsonSBDist:  possible lack of accuracy on variance");
      double C = -4.0 * Math.PI * delta * sum;
 
      double D = Math.sqrt(2 * Math.PI) * Math.exp(0.5 * gamma * gamma);
      double[] tpsi = new double[1];
      double mu = getMeanPsi(gamma, delta, xi, lambda, tpsi);
 
      double tvar = mu * (1 - delta * gamma) - mu * mu + delta / tpsi[0] * (A - B - mu * C * D);
      return lambda * lambda * tvar;
 
   }

   public static double getStandardDeviation(double gamma, double delta, double xi, double lambda) {
      return Math.sqrt(JohnsonSBDist.getVariance(gamma, delta, xi, lambda));
   }

   public void setParams(double gamma, double delta, double xi, double lambda) {
      setParams0(gamma, delta, xi, lambda);
      setLastParams(xi, lambda);
   }
 
}