Java Code Examples for org.apache.commons.math3.special.Gamma#digamma()

@Override
public void learn(final List<Datum> data, final double[] responsibilities) {
    double alpha = betaDistributionShape.getAlpha();
    double beta = betaDistributionShape.getBeta();

    for (int epoch = 0; epoch < NUM_EPOCHS; epoch++) {
        for (int n = 0; n < data.size(); n++) {
            final Datum datum = data.get(n);
            final int alt = datum.getAltCount();
            final int ref = datum.getTotalCount() - alt;

            final double digammaOfTotalPlusAlphaPlusBeta = Gamma.digamma(datum.getTotalCount() + alpha + beta);
            final double digammaOfAlphaPlusBeta = Gamma.digamma(alpha + beta);
            final double alphaGradient = Gamma.digamma(alpha + alt) - digammaOfTotalPlusAlphaPlusBeta - Gamma.digamma(alpha) + digammaOfAlphaPlusBeta;
            final double betaGradient = Gamma.digamma(beta + ref) - digammaOfTotalPlusAlphaPlusBeta - Gamma.digamma(beta) + digammaOfAlphaPlusBeta;

            alpha = Math.max(alpha + RATE * alphaGradient * responsibilities[n], 1.0);
            beta = Math.max(beta + RATE * betaGradient * responsibilities[n], 0.5);
        }
    }

    betaDistributionShape = new BetaDistributionShape(alpha, beta);

}
 

public void precomputeDenominatorForVariationalBayes( final double sumHyperParameterLambda ) {

        // Variational Bayes calculations from Bishop
        precomputeInverse();
        cachedSigmaInverse.timesEquals( hyperParameter_a );
        double sum = 0.0;
        for(int jjj = 1; jjj <= mu.length; jjj++) {
            sum += Gamma.digamma( (hyperParameter_a + 1.0 - jjj) / 2.0 );
        }
        sum -= Math.log( sigma.det() );
        sum += Math.log(2.0) * mu.length;
        final double lambda = 0.5 * sum;
        final double pi = Gamma.digamma( hyperParameter_lambda ) - Gamma.digamma( sumHyperParameterLambda );
        final double beta = (-1.0 * mu.length) / (2.0 * hyperParameter_b);
        cachedDenomLog10 = (pi / Math.log(10.0)) + (lambda / Math.log(10.0)) + (beta / Math.log(10.0));
    }
 

protected double[] calcEffectivePhis(final double E_alpha, final double[] responsibilitiesByRho) {

        final double sumResponsibilities = MathUtils.sum(responsibilitiesByRho);

        final double[] result = new double[responsibilitiesByRho.length];
        final int k = responsibilitiesByRho.length;

        // Initialize all pseudocounts to 1, except for index 0, which is 20;
        //  This artificially increases the odds for a rho = 0.
        final RealVector pseudocounts = new ArrayRealVector(responsibilitiesByRho.length);
        pseudocounts.set(1.0);
        pseudocounts.setEntry(0, 20.0);

        final double sumPseudocounts = MathUtils.sum(pseudocounts.toArray());
        final double term2 = Gamma.digamma(E_alpha + sumPseudocounts + sumResponsibilities);
        for (int i=0; i < result.length; i++) {
            final double term1 = Gamma.digamma(E_alpha/k + pseudocounts.getEntry(i) + responsibilitiesByRho[i]);
            result[i] = Math.exp(term1 - term2);
        }

        return result;
    }
 

@Override
public double[] calcExpectationLog() {
    double[] exLogMixing = new double[K];
    for (int i = 0; i < K; i++) {
        exLogMixing[i] = Gamma.digamma(coeffs[i]) - Gamma.digamma(sumCoeffs);
    }
    return exLogMixing;
}
 

/**
 * Overloaded digamma function: returns 0 if argument is zero
 *
 * @param x
 * @return digamma(x) if x nonzero, otherwise zero
 */
private double digamma(final double x) {
	if (Math.abs(x) < 1e-15)
		return 0.0;
	else
		return Gamma.digamma(x);
}
 

@Override
public double[] calcExpectationLog() {
    double[] lnMixingContribution = new double[T];
    double cumulativeAlreadyAssigned = 0;
    for (int t = 0; t < T; t++) {
        // Calculate Mixing coefficient contributions to r
        lnMixingContribution[t] = cumulativeAlreadyAssigned + (Gamma.digamma(shapeParams[t][0]) - Gamma.digamma(shapeParams[t][0] + shapeParams[t][1]));
        cumulativeAlreadyAssigned += Gamma.digamma(shapeParams[t][1]) - Gamma.digamma(shapeParams[t][0] + shapeParams[t][1]);
    }
    return lnMixingContribution;
}
 

public double getExpectationLogDeterminantLambda() {
    double ex_log_lambda = D * Math.log(2) + logDetOmega;
    for (int i=0; i<D; i++) {
        ex_log_lambda += Gamma.digamma(0.5 * (nu - i));
    }
    return ex_log_lambda;
}
 

private double expectationLnLambda() {
    double ex_ln_lambda = D * Math.log(2) + logDetOmega;
    for (int i = 1; i <= D; i++) {
        ex_ln_lambda += Gamma.digamma(0.5 * (nu + 1 - i));
    }
    return ex_ln_lambda;
}
 

private void updatePhiPerDoc(@Nonnegative final int d,
        @Nonnull final Map<String, float[]> eLogBeta_d) {
    // Dirichlet expectation (2d) for gamma
    final float[] gamma_d = _gamma[d];
    final double digamma_gammaSum_d = Gamma.digamma(MathUtils.sum(gamma_d));
    final double[] eLogTheta_d = new double[_K];
    for (int k = 0; k < _K; k++) {
        eLogTheta_d[k] = Gamma.digamma(gamma_d[k]) - digamma_gammaSum_d;
    }

    // updating phi w/ normalization
    final Map<String, float[]> phi_d = _phi.get(d);
    final Map<String, Float> doc = _miniBatchDocs.get(d);
    for (String label : doc.keySet()) {
        final float[] phi_label = phi_d.get(label);
        final float[] eLogBeta_label = eLogBeta_d.get(label);

        double normalizer = 0.d;
        for (int k = 0; k < _K; k++) {
            float phiVal = (float) Math.exp(eLogBeta_label[k] + eLogTheta_d[k]) + 1E-20f;
            phi_label[k] = phiVal;
            normalizer += phiVal;
        }

        for (int k = 0; k < _K; k++) {
            phi_label[k] /= normalizer;
        }
    }
}
 

@Nonnull
public static double[] digamma(@Nonnull final double[] arr) {
    final int k = arr.length;
    final double[] ret = new double[k];
    for (int i = 0; i < k; i++) {
        ret[i] = Gamma.digamma(arr[i]);
    }
    return ret;
}
 

@Nonnull
public static float[] digamma(@Nonnull final float[] arr) {
    final int k = arr.length;
    final float[] ret = new float[k];
    for (int i = 0; i < k; i++) {
        ret[i] = (float) Gamma.digamma(arr[i]);
    }
    return ret;
}
 

/**
 * Computes the inverse of the PSI function. Source:
 * http://bariskurt.com/calculating-the-inverse-of-digamma-function/
 *
 * @param value value
 * @return PSI^-1(value)
 */
private double inversePsi(double value) {
  double inv = (value >= -2.22) ? Math.exp(value) + 0.5 : -1.0 / (value - Gamma.digamma(1));

  for (int i = 0; i < 3; i++) {
    inv -= (Gamma.digamma(inv) - value) / Gamma.trigamma(inv);
  }

  return inv;
}
 

/**
 * Calculate digamma of the input value.
 */
public static double digamma(double x) {
    return Gamma.digamma(x);
}
 

/**
 * {@inheritDoc}
 */
@Override
public SufficientStatistics createInitSufficientStatistics() {
    return new SparseVectorDefaultValue(nOfStates,Gamma.digamma(1.0) - Gamma.digamma(this.sizeOfSufficientStatistics()));
}
 

@Override
protected double compute(final int n) {
    return Gamma.digamma(n);
}
 

public double[] effectiveMultinomialWeights() {
    final double digammaOfSum = Gamma.digamma(MathUtils.sum(alpha));
    return MathUtils.applyToArray(alpha, a -> Math.exp(Gamma.digamma(a) - digammaOfSum));
}
 

public double[] effectiveLog10MultinomialWeights() {
    final double digammaOfSum = Gamma.digamma(MathUtils.sum(alpha));
    return MathUtils.applyToArray(alpha, a -> (Gamma.digamma(a) - digammaOfSum) * MathUtils.LOG10_E);
}
 

public double[] effectiveLogMultinomialWeights() {
    final double digammaOfSum = Gamma.digamma(MathUtils.sum(alpha));
    return MathUtils.applyToArray(alpha, a -> (Gamma.digamma(a) - digammaOfSum));
}
 

@Override
public void fit() {
  System.out.println("\nFitting " + this.toString());

  for (int iter = 1; iter <= this.numIters; iter++) {
    Parallelizer.exec(this.datamodel.getUsers(), new UpdatePhiGamma());
    Parallelizer.exec(this.datamodel.getItems(), new UpdateBeta());

    double diff;

    do {
      diff = 0;

      double as = 0; // alpha sum
      for (int z = 0; z < this.numFactors; z++) {
        as += this.alpha[z];
      }

      for (int z = 0; z < this.numFactors; z++) {

        double gs = 0; // gamma sum
        for (int userIndex = 0; userIndex < datamodel.getNumberOfUsers(); userIndex++) {
          gs += this.gamma[userIndex][z];
        }

        double sum = 0;
        for (int userIndex = 0; userIndex < datamodel.getNumberOfUsers(); userIndex++) {
          sum += Gamma.digamma(this.gamma[userIndex][z]) - Gamma.digamma(gs);
        }

        double psiAlpha = Gamma.digamma(as) + sum / datamodel.getNumberOfUsers();

        double newAlpha = this.inversePsi(psiAlpha);

        diff += Math.abs(this.alpha[z] - newAlpha);

        this.alpha[z] = newAlpha;
      }
    } while (diff > EPSILON);

    if ((iter % 10) == 0) System.out.print(".");
    if ((iter % 100) == 0) System.out.println(iter + " iterations");
  }
}