Java FastMath.exp方法代码示例

import net.jafama.FastMath; //导入方法依赖的package包/类
/**
 * Convert Hoeffding D value to a p-value.
 * 
 * @param d D value
 * @param n Data set size
 * @return p-value
 */
public double toPValue(double d, int n) {
  double b = d / 30 + 1. / (36 * n);
  double z = .5 * MathUtil.PISQUARE * MathUtil.PISQUARE * n * b;

  // Exponential approximation
  if(z < 1.1 || z > 8.5) {
    double e = FastMath.exp(0.3885037 - 1.164879 * z);
    return (e > 1) ? 1 : (e < 0) ? 0 : e;
  }
  // Tabular approximation
  for(int i = 0; i < 86; i++) {
    if(TABPOS[i] >= z) {
      // Exact table value
      if(TABPOS[i] == z) {
        return TABVAL[i];
      }
      // Linear interpolation
      double x1 = TABPOS[i], x0 = TABPOS[i - 1];
      double y1 = TABVAL[i], y0 = TABVAL[i - 1];
      return y0 + (y1 - y0) * (z - x0) / (x1 - x0);
    }
  }
  return -1;
}
 

import net.jafama.FastMath; //导入方法依赖的package包/类
/**
 * Assigns the current probability values to the instances in the database and
 * compute the expectation value of the current mixture of distributions.
 * 
 * Computed as the sum of the logarithms of the prior probability of each
 * instance.
 * 
 * @param relation the database used for assignment to instances
 * @param models Cluster models
 * @param probClusterIGivenX Output storage for cluster probabilities
 * @return the expectation value of the current mixture of distributions
 */
public static double assignProbabilitiesToInstances(Relation<? extends NumberVector> relation, List<? extends EMClusterModel<?>> models, WritableDataStore<double[]> probClusterIGivenX) {
  final int k = models.size();
  double emSum = 0.;

  for(DBIDIter iditer = relation.iterDBIDs(); iditer.valid(); iditer.advance()) {
    NumberVector vec = relation.get(iditer);
    double[] probs = new double[k];
    for(int i = 0; i < k; i++) {
      probs[i] = models.get(i).estimateLogDensity(vec);
    }
    double logP = logSumExp(probs);
    emSum += logP > MIN_LOGLIKELIHOOD ? logP : MIN_LOGLIKELIHOOD;
    for(int i = 0; i < k; i++) {
      probs[i] = FastMath.exp(probs[i] - logP);
    }
    probClusterIGivenX.put(iditer, probs);
  }
  return emSum / relation.size();
}
 

import net.jafama.FastMath; //导入方法依赖的package包/类
/**
 * Compute H (observed perplexity) for row i, and the row pij_i.
 * 
 * @param ignore Object to skip
 * @param it Distances list
 * @param p Output probabilities
 * @param mbeta {@code -1. / (2 * sigma * sigma)}
 * @return Observed perplexity
 */
protected static double computeH(DBIDRef ignore, DoubleDBIDListIter it, double[] p, double mbeta) {
  double sumP = 0.;
  // Skip point "i", break loop in two:
  it.seek(0);
  for(int j = 0; it.valid(); j++, it.advance()) {
    if(DBIDUtil.equal(ignore, it)) {
      p[j] = 0;
      continue;
    }
    sumP += (p[j] = FastMath.exp(it.doubleValue() * mbeta));
  }
  if(!(sumP > 0)) {
    // All pij are zero. Bad news.
    return Double.NEGATIVE_INFINITY;
  }
  final double s = 1. / sumP; // Scaling factor
  double sum = 0.;
  // While we could skip pi[i], it should be 0 anyway.
  it.seek(0);
  for(int j = 0; it.valid(); j++, it.advance()) {
    sum += it.doubleValue() * (p[j] *= s);
  }
  return FastMath.log(sumP) - mbeta * sum;
}
 

import net.jafama.FastMath; //导入方法依赖的package包/类
@Override
public SkewGeneralizedNormalDistribution estimateFromLMoments(double[] xmom) {
  if(!(xmom[1] > 0.) || !(Math.abs(xmom[2]) < 1.0)) {
    throw new ArithmeticException("L-Moments invalid");
  }
  // Generalized Normal Distribution estimation:
  double t3 = xmom[2];
  if(Math.abs(t3) >= .95) {
    // Extreme skewness
    return new SkewGeneralizedNormalDistribution(0., -1., 0.);
  }
  if(Math.abs(t3) <= 1e-8) {
    // t3 effectively zero, this is a normal distribution.
    return new SkewGeneralizedNormalDistribution(xmom[0], xmom[1] * MathUtil.SQRTPI, 0.);
  }
  final double tt = t3 * t3;
  double shape = -t3 * (A0 + tt * (A1 + tt * (A2 + tt * A3))) / (1. + tt * (B1 + tt * (B2 + tt * B3)));
  final double e = FastMath.exp(.5 * shape * shape);
  double scale = xmom[1] * shape / (e * NormalDistribution.erf(.5 * shape));
  double location = xmom[0] + scale * (e - 1.) / shape;
  return new SkewGeneralizedNormalDistribution(location, scale, shape);
}
 

import net.jafama.FastMath; //导入方法依赖的package包/类
/**
 * Cumulative density function.
 * 
 * @param val Value
 * @param loc Location
 * @param scale Scale
 * @param shape Shape
 * @return CDF
 */
public static double cdf(double val, double loc, double scale, double shape) {
  val = (val - loc) / scale;
  if(shape != 0.) {
    final double tmp = 1 - shape * val;
    if(tmp < 1e-15) {
      return (shape < 0) ? 0 : 1;
    }
    val = -FastMath.log(tmp) / shape;
  }
  return 1. / (1. + FastMath.exp(-val));
}
 

import net.jafama.FastMath; //导入方法依赖的package包/类
@Override
public double pdf(double val) {
  if(val < location) {
    return 0.;
  }
  return rate * FastMath.exp(-rate * (val - location));
}
 

import net.jafama.FastMath; //导入方法依赖的package包/类
@Override
public double cdf(double val) {
  if(val < location) {
    return 0.;
  }
  return 1 - FastMath.exp(-rate * (val - location));
}
 

import net.jafama.FastMath; //导入方法依赖的package包/类
/**
 * Refinement of ChiSquared probit using Newton iterations.
 * 
 * A trick used by GNU R to improve precision.
 * 
 * @param logpt Target value of log p
 * @param k Alpha
 * @param theta Theta = 1 / Beta
 * @param maxit Maximum number of iterations to do
 * @param x Initial estimate
 * @return Refined value
 */
protected static double gammaQuantileNewtonRefinement(final double logpt, final double k, final double theta, final int maxit, double x) {
  final double EPS_N = 1e-15; // Precision threshold
  // 0 is not possible, try MIN_NORMAL instead
  if(x <= 0) {
    x = Double.MIN_NORMAL;
  }
  // Current estimation
  double logpc = logcdf(x, k, theta);
  if(x == Double.MIN_NORMAL && logpc > logpt * (1. + 1e-7)) {
    return 0.;
  }
  if(logpc == Double.NEGATIVE_INFINITY) {
    return 0.;
  }
  // Refine by newton iterations
  for(int i = 0; i < maxit; i++) {
    // Error of current approximation
    final double logpe = logpc - logpt;
    if(Math.abs(logpe) < Math.abs(EPS_N * logpt)) {
      break;
    }
    // Step size is controlled by PDF:
    final double g = logpdf(x, k, theta);
    if(g == Double.NEGATIVE_INFINITY) {
      break;
    }
    final double newx = x - logpe * FastMath.exp(logpc - g);
    // New estimate:
    logpc = logcdf(newx, k, theta);
    if(Math.abs(logpc - logpt) > Math.abs(logpe) || (i > 0 && Math.abs(logpc - logpt) == Math.abs(logpe))) {
      // no further improvement
      break;
    }
    x = newx;
  }
  return x;
}
 

import net.jafama.FastMath; //导入方法依赖的package包/类
/**
 * Get Gaussian Weight using standard deviation for scaling. max is ignored.
 */
@Override
public double getWeight(double distance, double max, double stddev) {
  if(stddev <= 0) {
    return 1;
  }
  double normdistance = distance / stddev;
  return scaling * FastMath.exp(-.5 * normdistance * normdistance) / stddev;
}
 

import net.jafama.FastMath; //导入方法依赖的package包/类
@Override
public GeneralizedExtremeValueDistribution estimateFromLMoments(double[] xmom) {
  double t3 = xmom[2];
  if(Math.abs(t3) < 1e-50 || (t3 >= 1.)) {
    throw new ArithmeticException("Invalid moment estimation.");
  }
  // Approximation for t3 between 0 and 1:
  double g;
  if(t3 > 0.) {
    double z = 1. - t3;
    g = (-1. + z * (C1 + z * (C2 + z * C3))) / (1. + z * (D1 + z * D2));
    // g: Almost zero?
    if(Math.abs(g) < 1e-50) {
      double k = 0;
      double sigma = xmom[1] * MathUtil.ONE_BY_LOG2;
      double mu = xmom[0] - EU * sigma;
      return new GeneralizedExtremeValueDistribution(mu, sigma, k);
    }
  }
  else {
    // Approximation for t3 between -.8 and 0L:
    g = (A0 + t3 * (A1 + t3 * (A2 + t3 * (A3 + t3 * A4)))) / (1. + t3 * (B1 + t3 * (B2 + t3 * B3)));
    if(t3 < -.8) {
      // Newton-Raphson iteration for t3 < -.8
      if(t3 <= -.97) {
        g = 1. - FastMath.log1p(t3) * MathUtil.ONE_BY_LOG2;
      }
      double t0 = .5 * (t3 + 3.);
      for(int it = 1;; it++) {
        // These sometimes do not converge with FastMath.
        double x2 = FastMath.pow(2., -g), xx2 = 1. - x2;
        double x3 = FastMath.pow(3., -g), xx3 = 1. - x3;
        double t = xx3 / xx2;
        double deriv = (xx2 * x3 * MathUtil.LOG3 - xx3 * x2 * MathUtil.LOG2) / (xx2 * xx2);
        double oldg = g;
        g -= (t - t0) / deriv;
        if(Math.abs(g - oldg) < 1e-14 * g) {
          break;
        }
        if(it >= MAXIT) {
          throw new ArithmeticException("Newton-Raphson did not converge.");
        }
      }
    }
  }
  double gam = FastMath.exp(GammaDistribution.logGamma(1. + g));
  final double mu, sigma, k;
  k = g;
  sigma = xmom[1] * g / (gam * (1. - FastMath.pow(2., -g)));
  mu = xmom[0] - sigma * (1. - gam) / g;
  return new GeneralizedExtremeValueDistribution(mu, sigma, k);
}
 

import net.jafama.FastMath; //导入方法依赖的package包/类
/**
 * Run the algorithm
 * 
 * @param relation Data relation
 * @return Outlier result
 */
public OutlierResult run(Relation<V> relation) {
  DoubleMinMax mm = new DoubleMinMax();
  // resulting scores
  WritableDoubleDataStore oscores = DataStoreUtil.makeDoubleStorage(relation.getDBIDs(), DataStoreFactory.HINT_TEMP | DataStoreFactory.HINT_HOT);

  // Compute mean and covariance Matrix
  CovarianceMatrix temp = CovarianceMatrix.make(relation);
  double[] mean = temp.getMeanVector(relation).toArray();
  // debugFine(mean.toString());
  double[][] covarianceMatrix = temp.destroyToPopulationMatrix();
  // debugFine(covarianceMatrix.toString());
  double[][] covarianceTransposed = inverse(covarianceMatrix);

  // Normalization factors for Gaussian PDF
  double det = new LUDecomposition(covarianceMatrix).det();
  final double fakt = 1.0 / FastMath.sqrt(MathUtil.powi(MathUtil.TWOPI, RelationUtil.dimensionality(relation)) * det);

  // for each object compute Mahalanobis distance
  for(DBIDIter iditer = relation.iterDBIDs(); iditer.valid(); iditer.advance()) {
    double[] x = minusEquals(relation.get(iditer).toArray(), mean);
    // Gaussian PDF
    final double mDist = transposeTimesTimes(x, covarianceTransposed, x);
    final double prob = fakt * FastMath.exp(-mDist * .5);

    mm.put(prob);
    oscores.putDouble(iditer, prob);
  }

  final OutlierScoreMeta meta;
  if(invert) {
    double max = mm.getMax() != 0 ? mm.getMax() : 1.;
    for(DBIDIter iditer = relation.iterDBIDs(); iditer.valid(); iditer.advance()) {
      oscores.putDouble(iditer, (max - oscores.doubleValue(iditer)) / max);
    }
    meta = new BasicOutlierScoreMeta(0.0, 1.0);
  }
  else {
    meta = new InvertedOutlierScoreMeta(mm.getMin(), mm.getMax(), 0.0, Double.POSITIVE_INFINITY);
  }
  DoubleRelation res = new MaterializedDoubleRelation("Gaussian Model Outlier Score", "gaussian-model-outlier", oscores, relation.getDBIDs());
  return new OutlierResult(meta, res);
}
 

import net.jafama.FastMath; //导入方法依赖的package包/类
/**
 * Complementary error function for Gaussian distributions = Normal
 * distributions.
 * 
 * Based on:<br />
 * Takuya Ooura, http://www.kurims.kyoto-u.ac.jp/~ooura/gamerf.html<br />
 * Copyright (C) 1996 Takuya OOURA (email: [email protected]).<br />
 * "You may use, copy, modify this code for any purpose and without fee."
 * 
 * @param x parameter value
 * @return erfc(x)
 */
@Reference(authors = "T. Ooura", //
    title = "Gamma / Error Functions", //
    booktitle = "Online", //
    url = "http://www.kurims.kyoto-u.ac.jp/~ooura/gamerf.html")
public static double erfc(double x) {
  if(Double.isNaN(x)) {
    return Double.NaN;
  }
  if(Double.isInfinite(x)) {
    return (x < 0.0) ? 2 : 0;
  }
  final double t = 3.97886080735226 / (Math.abs(x) + 3.97886080735226);
  final double u = t - 0.5;
  double y = (((//
  ((((((0.00127109764952614092 * u //
      + 1.19314022838340944e-4) * u //
      - 0.003963850973605135) * u //
      - 8.70779635317295828e-4) * u //
      + 0.00773672528313526668) * u //
      + 0.00383335126264887303) * u //
      - 0.0127223813782122755) * u //
      - 0.0133823644533460069) * u //
      + 0.0161315329733252248) * u //
      + 0.0390976845588484035) * u //
      + 0.00249367200053503304;
  y = ((((((((((((y * u //
      - 0.0838864557023001992) * u //
      - 0.119463959964325415) * u //
      + 0.0166207924969367356) * u //
      + 0.357524274449531043) * u //
      + 0.805276408752910567) * u //
      + 1.18902982909273333) * u //
      + 1.37040217682338167) * u //
      + 1.31314653831023098) * u //
      + 1.07925515155856677) * u //
      + 0.774368199119538609) * u //
      + 0.490165080585318424) * u //
      + 0.275374741597376782) //
      * t * FastMath.exp(-x * x);
  return x < 0 ? 2 - y : y;
}
 

import net.jafama.FastMath; //导入方法依赖的package包/类
@Override
public double pdf(double val) {
  return .5 * rate * FastMath.exp(-rate * Math.abs(val - location));
}
 

import net.jafama.FastMath; //导入方法依赖的package包/类
/**
 * PDF of Gumbel distribution
 * 
 * @param x Value
 * @param mu Mode
 * @param beta Shape
 * @return PDF at position x.
 */
public static double pdf(double x, double mu, double beta) {
  final double z = (x - mu) / beta;
  if(x == Double.NEGATIVE_INFINITY) {
    return 0.;
  }
  return FastMath.exp(-z - FastMath.exp(-z)) / beta;
}
 

import net.jafama.FastMath; //导入方法依赖的package包/类
/**
 * Probability density function of the ExGaussian distribution.
 *
 * @param x The value.
 * @param mu The mean.
 * @param sigma The standard deviation.
 * @param lambda Rate parameter.
 * @return PDF of the given exgauss distribution at x.
 */
public static double pdf(double x, double mu, double sigma, double lambda) {
  final double dx = x - mu;
  final double lss = lambda * sigma * sigma;
  final double erfc = NormalDistribution.erfc((lss - dx) / (sigma * MathUtil.SQRT2));
  return erfc > 0. ? .5 * lambda * FastMath.exp(lambda * (lss * .5 - dx)) * erfc : 0.;
}