Example usage for org.apache.commons.math.special Gamma logGamma

List of usage examples for org.apache.commons.math.special Gamma logGamma

Introduction

In this page you can find the example usage for org.apache.commons.math.special Gamma logGamma.

Prototype

public static double logGamma(double x) 

Source Link

Document

Returns the natural logarithm of the gamma function Γ(x).

Usage

From source file:com.opengamma.analytics.math.function.special.NaturalLogGammaFunction.java

/**
 * @param x The argument of the function, must be greater than zero
 * @return The value of the function // w  w  w .  ja  v a 2s .c o  m
 */
@Override
public Double evaluate(final Double x) {
    Validate.isTrue(x > 0, "x must be greater than zero");
    return Gamma.logGamma(x);
}

From source file:com.opengamma.analytics.math.function.special.GammaFunction.java

@Override
public Double evaluate(final Double x) {
    if (x > 0.0) {
        return Math.exp(Gamma.logGamma(x));
    }/*from  ww  w  .j ava  2 s  .  c o  m*/
    return Math.PI / Math.sin(Math.PI * x) / evaluate(1 - x);
}

From source file:com.cloudera.science.mgps.NFunction.java

public double eval(int n, double e) {
    double x = -n * Math.log(1 + beta / e);
    double y = -alpha * Math.log(1 + e / beta);
    double z = Gamma.logGamma(alpha + n);
    double d = Gamma.logGamma(alpha) + MathUtils.factorialLog(n);
    return Math.exp(x + y + z - d);
}

From source file:com.analog.lyric.chimple.monkeys.ChimpDirichlet.java

@Override
public double calculateLogLikelihood(Object result, Object[] parameters) {
    double[] alphas = (double[]) parameters[0];
    double[] value = (double[]) result;

    double retval = 0;

    double sum = 0;
    for (int i = 0; i < alphas.length; i++) {
        retval -= (alphas[i] - 1) * Math.log(value[i]);
        retval += Gamma.logGamma(alphas[i]);
        sum += alphas[i];//w w  w. j a va2 s.c  o m
    }

    retval -= Gamma.logGamma(sum);

    //result=sum(gammaln(alphas))-gammaln(sum(alphas))-sum((alphas-1).*log(value));

    return retval;
}

From source file:com.analog.lyric.chimple.monkeys.ChimpBeta.java

@Override
public double calculateLogLikelihood(Object result, Object[] parameters) {
    /*/*from  w  w  w .ja  va  2 s. c om*/
     *     alphas=[alpha beta];
    value=[value 1-value];
    %fast computation without normalization
    %result=-sum((alphas-1).*log(value));
            
    %exact computation, but the constant term should hopefully be useless
            
    result=sum(gammaln(alphas))-gammaln(sum(alphas))-sum((alphas-1).*log(value));
            
     */
    double value = (Double) result;
    double alpha = (Double) parameters[0];
    double beta = (Double) parameters[1];

    //double alphaBetaSum = alpha+beta;

    double retval = Gamma.logGamma(alpha);
    retval += Gamma.logGamma(beta);
    retval -= Gamma.logGamma(alpha + beta);
    retval -= (alpha - 1) * Math.log(value);
    retval -= (beta - 1) * Math.log(1 - value);

    return retval;
}

From source file:geogebra.util.MyMath.java

final public static double gamma(double x, Kernel kernel) {

    // Michael Borcherds 2008-05-04
    if (x <= 0 && Kernel.isEqual(x, Math.round(x)))
        return Double.NaN; // negative integers

    // Michael Borcherds 2007-10-15 BEGIN added case for x<0 otherwise no
    // results in 3rd quadrant
    if (x >= 0)
        return Math.exp(Gamma.logGamma(x));
    else//from  w w w. j  av  a 2  s.c  o m
        return -Math.PI / (x * gamma(-x, kernel) * Math.sin(Math.PI * x));
    // Michael Borcherds 2007-10-15 END
}

From source file:com.opengamma.analytics.math.statistics.distribution.NonCentralChiSquaredDistribution.java

/**
 * @param degrees The number of degrees of freedom, not negative or zero
 * @param nonCentrality The non-centrality parameter, not negative
 *///w ww.  j a  v  a  2s  .  c  o  m
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);
    }
}

From source file:master.steppers.TauLeapingStepper.java

/**
 * Generate appropriate random state change according to Gillespie's
 * tau-leaping algorithm.//  www .ja  v a 2s. co m
 *
 * @param reaction
 * @param state PopulationState to modify.
 * @param model
 * @param calcLogP
 * @param thisdt
 */
public void leap(Reaction reaction, PopulationState state, Model model, boolean calcLogP, double thisdt) {

    // Draw number of reactions to fire within time tau:
    double rho = reaction.getPropensity() * thisdt;
    double q = Randomizer.nextPoisson(rho);

    if (calcLogP) {
        if (rho > 0)
            stepLogP += -rho + q * Math.log(rho / thisdt) - Gamma.logGamma(q + 1);
    }

    // Implement reactions:
    state.implementReaction(reaction, q);

    // Increment event counter:
    eventCount += q;

}

From source file:dr.math.distributions.BetaDistribution.java

private void recomputeZ() {
    if (Double.isNaN(z)) {
        z = Gamma.logGamma(alpha) + Gamma.logGamma(beta) - Gamma.logGamma(alpha + beta);
    }
}

From source file:master.steppers.SALStepper.java

/**
 * Generate appropriate random state change according to Sehl's
 * step anticipation tau-leaping algorithm.
 *
 * @param reaction/* www . j  a v a2s.co m*/
 * @param state PopulationState to modify.
 * @param model
 * @param calcLogP
 * @param thisdt
 */
public void leap(Reaction reaction, PopulationState state, Model model, boolean calcLogP, double thisdt) {

    // Calculate corrected rate
    double rho = reaction.getPropensity() * thisdt + 0.5 * corrections.get(reaction) * thisdt * thisdt;

    // Draw number of reactions to fire within time tau:
    double q = Randomizer.nextPoisson(rho);

    if (calcLogP) {
        if (rho > 0)
            stepLogP += -rho + q * Math.log(rho / thisdt) - Gamma.logGamma(q + 1);
    }

    // Implement reactions:
    state.implementReaction(reaction, q);

    // Increment event counter:
    eventCount += q;

}