InverseIncompleteBetaFunction.java

/**
 * Copyright (C) 2009 - present by OpenGamma Inc. and the OpenGamma group of companies
 *
 * Please see distribution for license.
 */
package com.opengamma.strata.math.impl.function.special;

import java.util.function.Function;

import com.google.common.math.DoubleMath;
import com.opengamma.strata.collect.ArgChecker;
import com.opengamma.strata.math.MathException;

/**
 * 
 */
public class InverseIncompleteBetaFunction implements Function<Double, Double> {
//TODO either find another implementation or delete this class

  private final double _a;
  private final double _b;
  private final Function<Double, Double> _lnGamma = new NaturalLogGammaFunction();
  private final Function<Double, Double> _beta;
  private static final double EPS = 1e-9;

  public InverseIncompleteBetaFunction(double a, double b) {
    ArgChecker.notNegativeOrZero(a, "a");
    ArgChecker.notNegativeOrZero(b, "b");
    _a = a;
    _b = b;
    _beta = new IncompleteBetaFunction(a, b);
  }

  //-------------------------------------------------------------------------
  @Override
  public Double apply(Double x) {
    ArgChecker.inRangeInclusive(x, 0d, 1d, "x");
    double pp, p, t, h, w, lnA, lnB, u, a1 = _a - 1;
    double b1 = _b - 1;
    if (_a >= 1 && _b >= 1) {
      pp = x < 0.5 ? x : 1 - x;
      t = Math.sqrt(-2 * Math.log(pp));
      p = (2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t;
      if (p < 0.5) {
        p *= -1;
      }
      a1 = (Math.sqrt(p) - 3.) / 6.;
      double tempA = 1. / (2 * _a - 1);
      double tempB = 1. / (2 * _b - 1);
      h = 2. / (tempA + tempB);
      w = p * Math.sqrt(a1 + h) / h - (tempB - tempA) * (a1 + 5. / 6 - 2. / (3 * h));
      p = _a / (_a + _b + Math.exp(2 * w));
    } else {
      lnA = Math.log(_a / (_a + _b));
      lnB = Math.log(_b / (_a + _b));
      t = Math.exp(_a * lnA) / _a;
      u = Math.exp(_b * lnB) / _b;
      w = t + u;
      if (x < t / w) {
        p = Math.pow(_a * w * x, 1. / _a);
      } else {
        p = 1 - Math.pow(_b * w * (1 - x), 1. / _b);
      }
    }
    double afac = -_lnGamma.apply(_a) - _lnGamma.apply(_b) + _lnGamma.apply(_a + _b);
    double error;
    for (int j = 0; j < 10; j++) {
      if (DoubleMath.fuzzyEquals(p, 0d, 1e-16) || DoubleMath.fuzzyEquals(p, (double) 1, 1e-16)) {
        throw new MathException("a or b too small for accurate evaluation");
      }
      error = _beta.apply(p) - x;
      t = Math.exp(a1 * Math.log(p) + b1 * Math.log(1 - p) + afac);
      u = error / t;
      t = u / (1 - 0.5 * Math.min(1, u * (a1 / p - b1 / (1 - p))));
      p -= t;
      if (p <= 0) {
        p = 0.5 * (p + t);
      }
      if (p >= 1) {
        p = 0.5 * (p + t + 1);
      }
      if (Math.abs(t) < EPS * p && j > 0) {
        break;
      }
    }
    return p;
  }

}