code: public static double quantile(double p, double k, double theta, double shift) {
      return FastMath.log(GammaDistribution.quantile(p, k, 1)) / theta + shift;
    }
docs: Compute probit (inverse cdf) for ExpGamma distributions.
  
  @param p Probability
  @param k k, alpha aka. "shape" parameter
  @param theta Theta = 1.0/Beta aka. "scaling" parameter
  @param shift Shift parameter
  @return Probit for ExpGamma distribution