import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
@Override
public PartitionGetResult partitionGet(PartitionGetParam partParam) {
  PartitionKey pkey = partParam.getPartKey();

  pkey = psContext.getMatrixMetaManager().getMatrixMeta(pkey.getMatrixId())
    .getPartitionMeta(pkey.getPartitionId()).getPartitionKey();

  int ws = pkey.getStartRow();
  int es = pkey.getEndRow();

  LikelihoodParam.LikelihoodPartParam param = (LikelihoodParam.LikelihoodPartParam) partParam;
  float beta = param.getBeta();

  double lgammaBeta = Gamma.logGamma(beta);

  double ll = 0;
  for (int w = ws; w < es; w ++) {
    ServerRow row = psContext.getMatrixStorageManager().getRow(pkey, w);
    ll += likelihood(row, beta, lgammaBeta);
  }

  return new ScalarPartitionAggrResult(ll);
}
 

import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
public Sampler(CSRTokens data, LDAModel model) {
  this.data = data;
  this.model = model;

  K = model.K();
  alpha = model.alpha();
  beta  = model.beta();
  vbeta = data.n_words * beta;

  lgammaBeta = Gamma.logGamma(beta);
  lgammaAlpha = Gamma.logGamma(alpha);
  lgammaAlphaSum = Gamma.logGamma(alpha * K);

  nk = new int[K];
  wk = new int[K];
  tidx = new short[K];
  psum = new float[K];

  tree = new FTree(K);
}
 

import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
public static double dpois_raw(double x, double lambda, boolean give_log) {
  /*       x >= 0 ; integer for dpois(), but not e.g. for pgamma()!
  lambda >= 0
   */
  if (lambda == 0) {
    return ((x == 0) ? SignRank.R_D__1(true, give_log) : SignRank.R_D__0(true, give_log));
  }
  if (!DoubleVector.isFinite(lambda)) {
    return SignRank.R_D__0(true, give_log);
  }
  if (x < 0) {
    return (SignRank.R_D__0(true, give_log));
  }
  if (x <= lambda * Double.MIN_VALUE) {
    return (SignRank.R_D_exp(-lambda, true, give_log));
  }
  if (lambda < x * Double.MIN_VALUE) {
    return (SignRank.R_D_exp(-lambda + x * Math.log(lambda) - Gamma.logGamma(x + 1), true, give_log));
  }
  return (SignRank.R_D_fexp(2 * Math.PI * x, -Binom.stirlerr(x) - Binom.bd0(x, lambda), true, give_log));
}
 

import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
public void testMomonts() {
    final double tol = 1e-9;
    WeibullDistribution dist;
    
    dist = new WeibullDistributionImpl(2.5, 3.5);
    // In R: 3.5*gamma(1+(1/2.5)) (or emperically: mean(rweibull(10000, 2.5, 3.5)))
    assertEquals(dist.getNumericalMean(), 3.5 * FastMath.exp(Gamma.logGamma(1 + (1 / 2.5))), tol);
    assertEquals(dist.getNumericalVariance(), (3.5 * 3.5) * 
            FastMath.exp(Gamma.logGamma(1 + (2 / 2.5))) -
            (dist.getNumericalMean() * dist.getNumericalMean()), tol); 
    
    dist = new WeibullDistributionImpl(10.4, 2.222);
    assertEquals(dist.getNumericalMean(), 2.222 * FastMath.exp(Gamma.logGamma(1 + (1 / 10.4))), tol);
    assertEquals(dist.getNumericalVariance(), (2.222 * 2.222) * 
            FastMath.exp(Gamma.logGamma(1 + (2 / 10.4))) -
            (dist.getNumericalMean() * dist.getNumericalMean()), tol);
}
 

import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
@Test
public void testMomonts() {
    final double tol = 1e-9;
    WeibullDistribution dist;
    
    dist = new WeibullDistributionImpl(2.5, 3.5);
    // In R: 3.5*gamma(1+(1/2.5)) (or emperically: mean(rweibull(10000, 2.5, 3.5)))
    Assert.assertEquals(dist.getNumericalMean(), 3.5 * FastMath.exp(Gamma.logGamma(1 + (1 / 2.5))), tol);
    Assert.assertEquals(dist.getNumericalVariance(), (3.5 * 3.5) * 
            FastMath.exp(Gamma.logGamma(1 + (2 / 2.5))) -
            (dist.getNumericalMean() * dist.getNumericalMean()), tol); 
    
    dist = new WeibullDistributionImpl(10.4, 2.222);
    Assert.assertEquals(dist.getNumericalMean(), 2.222 * FastMath.exp(Gamma.logGamma(1 + (1 / 10.4))), tol);
    Assert.assertEquals(dist.getNumericalVariance(), (2.222 * 2.222) * 
            FastMath.exp(Gamma.logGamma(1 + (2 / 10.4))) -
            (dist.getNumericalMean() * dist.getNumericalMean()), tol);
}
 

import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
protected double functionLogZ(double[] vector){

	double valueFunction;
	double sum=0;
	double logpart1=0;

	// log(prod(X))=log(X1)+log(X2)+.........+log(Xn)	
	for(int i=0;i<vector.length;i++){
		logpart1=logpart1+Gamma.logGamma(vector[i]);	
		sum=sum+vector[i];
	}	

	double logpart2 = Gamma.logGamma(sum);

	valueFunction = logpart1-logpart2;

	return valueFunction;
}
 

import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
private double likelihood(ServerRow row, float beta, double lgammaBeta) {
  int len = (int)(row.getEndCol() - row.getStartCol());
  double ll = 0;
  if (row instanceof ServerDenseIntRow) {
    IntBuffer buf = ((ServerDenseIntRow) row).getData();
    for (int i = 0; i < len; i ++) {
      if (buf.get(i) > 0)
        ll += Gamma.logGamma(buf.get(i) + beta) - lgammaBeta;
    }
  } else
    throw new AngelException("should be ServerDenseIntRow");
  return ll;
}
 

import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
@Override
public Double evaluate(final Double x) {
  try {
    return Gamma.regularizedGammaP(_a, x, _eps, _maxIter);
  } catch (final org.apache.commons.math.MathException e) {
    throw new MathException(e);
  }
}
 

import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
@Override
public Double evaluate(final Double x) {
  if (x > 0.0) {
    return Math.exp(Gamma.logGamma(x));
  }
  return Math.PI / Math.sin(Math.PI * x) / evaluate(1 - x);
}
 

import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
/**
 * @param degrees The number of degrees of freedom, not negative or zero
 * @param nonCentrality The non-centrality parameter, not negative
 */
public NonCentralChiSquaredDistribution(final double degrees, final double nonCentrality) {
  Validate.isTrue(degrees > 0, "degrees of freedom must be > 0, have " + degrees);
  Validate.isTrue(nonCentrality >= 0, "non-centrality must be >= 0, have " + nonCentrality);
  _dofOverTwo = degrees / 2.0;
  _lambdaOverTwo = nonCentrality / 2.0;
  _k = (int) Math.round(_lambdaOverTwo);

  if (_lambdaOverTwo == 0) {
    _pStart = 0.0;
  } else {
    final double logP = -_lambdaOverTwo + _k * Math.log(_lambdaOverTwo) - Gamma.logGamma(_k + 1);
    _pStart = Math.exp(logP);
  }
}
 

import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
/**
 * The probability distribution function P(X <= x) for a Poisson distribution.
 * 
 * @param x the value at which the PDF is evaluated.
 * @return Poisson distribution function evaluated at x
 * @throws MathException if the cumulative probability can not be
 *            computed due to convergence or other numerical errors.
 */
public double cumulativeProbability(int x) throws MathException {
    if (x < 0) {
        return 0;
    }
    if (x == Integer.MAX_VALUE) {
        return 1;
    }
    return Gamma.regularizedGammaQ((double)x + 1, mean, 
            1E-12, Integer.MAX_VALUE);
}
 

import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
/**
 * Returns the probability density for a particular point.
 *
 * @param x The point at which the density should be computed.
 * @return The pdf at point x.
 * @since 2.1
 */
@Override
public double density(double x) {
    final double n = degreesOfFreedom;
    final double plus1Over2 = (n + 1) / 2;
    return Math.exp(Gamma.logGamma(plus1Over2) - 0.5 * (Math.log(Math.PI) + Math.log(n)) -
            Gamma.logGamma(n / 2) - plus1Over2 * Math.log(1 + x * x / n));
}
 

import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
/**
 * The probability distribution function P(X <= x) for a Poisson
 * distribution.
 *
 * @param x the value at which the PDF is evaluated.
 * @return Poisson distribution function evaluated at x
 * @throws MathException if the cumulative probability can not be computed
 *                       due to convergence or other numerical errors.
 */
@Override
public double cumulativeProbability(int x) throws MathException {
    if (x < 0) {
        return 0;
    }
    if (x == Integer.MAX_VALUE) {
        return 1;
    }
    return Gamma.regularizedGammaQ((double) x + 1, mean, epsilon, maxIterations);
}
 

import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
public static Matrix gammaLn(Matrix m){

        int nc = m.data.length;
        int nr = m.data[0].length;
        double[][] result = new double[nc][];
        for (int k = 0; k < nc; ++k) {
            result[k] = new double[nr];
            for (int w = 0; w < nr; ++w) {
                result[k][w] = Gamma.logGamma(m.data[k][w]);
            }
        }
        return new Matrix(result);
    }
 

import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
/**
 * Digamma function (the first derivative of the logarithm of the gamma function).
 * @param array   - variational parameter
 * @return
 */
static double[][] dirichletExpectation(double[][] array) {
    int numRows = array.length;
    int numCols = array[0].length;

    double[] vector = new double[numRows];
    Arrays.fill(vector, 0.0);

    for (int k = 0; k < numRows; ++k) {
        for (int w = 0; w < numCols; ++w) {
            try{
                vector[k] += array[k][w];
            }catch (Exception e){
                throw new RuntimeException(e);
            }
        }
    }
    for (int k = 0; k < numRows; ++k) {
        vector[k] = Gamma.digamma(vector[k]);
    }

    double [][] approx = new double[numRows][];
    for (int k = 0; k < numRows; ++k) {
        approx[k] = new double[numCols];
        for (int w = 0; w < numCols; ++w) {
            double z =  Gamma.digamma(array[k][w]);
            approx[k][w] = z - vector[k];
        }
    }
    return approx;
}
 

import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
/**
 * The probability distribution function P(X <= x) for a Poisson
 * distribution.
 *
 * @param x the value at which the PDF is evaluated.
 * @return Poisson distribution function evaluated at x
 * @throws MathException if the cumulative probability can not be computed
 *             due to convergence or other numerical errors.
 */
@Override
public double cumulativeProbability(int x) throws MathException {
    if (x < 0) {
        return 0;
    }
    if (x == Integer.MAX_VALUE) {
        return 1;
    }
    return Gamma.regularizedGammaQ((double) x + 1, mean, 1E-12,
            Integer.MAX_VALUE);
}
 

import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
/**
 * {@inheritDoc}
 *
 * The variance is
 * <code>scale^2 * Gamma(1 + (2 / shape)) - mean^2</code>
 * where <code>Gamma(...)</code> is the Gamma-function
 *
 * @return {@inheritDoc}
 */
@Override
protected double calculateNumericalVariance() {
    final double sh = getShape();
    final double sc = getScale();
    final double mn = getNumericalMean();

    return (sc * sc) *
        FastMath.exp(Gamma.logGamma(1 + (2 / sh))) -
        (mn * mn);
}
 

import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
/**
 * {@inheritDoc}
 */
@Override
public double density(double x) {
    final double n = degreesOfFreedom;
    final double nPlus1Over2 = (n + 1) / 2;
    return FastMath.exp(Gamma.logGamma(nPlus1Over2) -
                        0.5 * (FastMath.log(FastMath.PI) + FastMath.log(n)) -
                        Gamma.logGamma(n/2) - nPlus1Over2 * FastMath.log(1 + x * x /n));
}
 

import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
/**
 * {@inheritDoc}
 */
@Override
public double density(double x) {
    if (x < 0) return 0;
    return FastMath.pow(x / beta, alpha - 1) / beta *
        FastMath.exp(-x / beta) / FastMath.exp(Gamma.logGamma(alpha));
}
 

import org.apache.commons.math.special.Gamma; //导入依赖的package包/类
/**
 * The probability distribution function {@code P(X <= x)} for a Poisson
 * distribution.
 *
 * @param x Value at which the PDF is evaluated.
 * @return the Poisson distribution function evaluated at {@code x}.
 * @throws MathException if the cumulative probability cannot be computed
 * due to convergence or other numerical errors.
 */
@Override
public double cumulativeProbability(int x) throws MathException {
    if (x < 0) {
        return 0;
    }
    if (x == Integer.MAX_VALUE) {
        return 1;
    }
    return Gamma.regularizedGammaQ((double) x + 1, mean, epsilon, maxIterations);
}