List of usage examples for org.apache.commons.math3.special Gamma trigamma
public static double trigamma(double x)
From source file:gedi.util.math.stat.distributions.LfcDistribution.java
public static double var(double a, double b) { return (Gamma.trigamma(a) + Gamma.trigamma(b)) / Math.log(2) / Math.log(2); }
From source file:edu.byu.nlp.stats.DirichletMLEOptimizable.java
/** {@inheritDoc} */ @Override/* w ww. j a v a 2 s .c om*/ public ValueAndObject<double[]> computeNext(double[] alpha) { if (!inPlace) { alpha = alpha.clone(); } double alphaSum = DoubleArrays.sum(alpha); double[] g = computeG(alpha, alphaSum); double[] q = computeQ(alpha); double b = computeB(g, q, N * Gamma.trigamma(alphaSum)); for (int k = 0; k < alpha.length; k++) { alpha[k] -= (g[k] - b) / q[k]; } double value = computeValue(alpha); return new ValueAndObject<double[]>(value, alpha); }
From source file:edu.byu.nlp.stats.SymmetricDirichletMLEOptimizableTest.java
/** * Test method for {@link edu.byu.nlp.stats.SymmetricDirichletMLENROptimizable#computeNext(java.lang.Double)}. *///from w ww . j a v a2 s.c om @Test public void testComputeNext() { final double alpha = 2.13; final double[][] data = DirichletTestUtils.sampleDataset(); SymmetricDirichletMLENROptimizable o = SymmetricDirichletMLENROptimizable.newOptimizable(data); ValueAndObject<Double> vao = o.computeNext(alpha); // Compute the update that should have happened final int N = data.length; final int K = data[0].length; double sumOfLogThetas = DoubleArrays.sum(Matrices.sumOverFirst(data)); double value = N * Gamma.logGamma(K * alpha) - N * K * Gamma.logGamma(alpha) + (alpha - 1) * sumOfLogThetas; double firstDeriv = N * Gamma.digamma(K * alpha) - N * K * Gamma.digamma(alpha) + sumOfLogThetas; double secondDeriv = N * Gamma.trigamma(K * alpha) - N * K * Gamma.trigamma(alpha); double nextAlpha = alpha - firstDeriv / secondDeriv; // Ensure that the update was computed correctly assertThat(vao.getValue()).isEqualTo(value, delta); assertThat(vao.getObject().doubleValue()).isEqualTo(nextAlpha, delta); }
From source file:edu.byu.nlp.stats.DirichletMLEOptimizableTest.java
/** * Computes a Newton-Raphson update in-place to alpha. */// w ww . jav a2 s . c o m private static RealVector newtonRaphsonUpdate(final double[][] data, double[] alpha) { // We'll compute the gold-standard value the "long" way (taking the inverse of the Hessian) RealMatrix hessian = new Array2DRowRealMatrix(alpha.length, alpha.length); for (int r = 0; r < hessian.getRowDimension(); r++) { for (int c = 0; c < hessian.getColumnDimension(); c++) { hessian.addToEntry(r, c, data.length * Gamma.trigamma(DoubleArrays.sum(alpha))); if (r == c) { hessian.addToEntry(r, c, -data.length * Gamma.trigamma(alpha[r])); } } } RealVector derivative = new ArrayRealVector(alpha.length); for (int k = 0; k < alpha.length; k++) { derivative.setEntry(k, data.length * (Gamma.digamma(DoubleArrays.sum(alpha)) - Gamma.digamma(alpha[k]))); for (double[] theta : data) { derivative.addToEntry(k, theta[k]); } } RealMatrix hessianInverse = new LUDecomposition(hessian).getSolver().getInverse(); RealVector negDiff = hessianInverse.preMultiply(derivative); negDiff.mapMultiplyToSelf(-1.0); RealVector expected = new ArrayRealVector(alpha, true); return expected.add(negDiff); }
From source file:edu.byu.nlp.stats.SymmetricDirichletOptimizationHelper.java
public double secondDerivativeAt(double alpha) { Preconditions.checkArgument(alpha > 0.0, "alpha must be strictly greater than zero; was %s", alpha); double secondDerivative = N * K * (K * Gamma.trigamma(K * alpha) - Gamma.trigamma(alpha)); logger.debug(String.format("secondDerivativeAt(%f) = %f", alpha, secondDerivative)); return secondDerivative; }
From source file:edu.byu.nlp.stats.DirichletMLEOptimizable.java
private double[] computeQ(double[] alpha) { double[] q = new double[alpha.length]; for (int k = 0; k < q.length; k++) { q[k] = -N * Gamma.trigamma(alpha[k]); }//from w ww.j av a2s.c o m return q; }
From source file:gamlss.distributions.GA.java
/** Second derivative d2lds2= (d^2l)/(dsigma^2), * where l - log-likelihood function./*from w ww . j a v a 2 s .c o m*/ * @return a vector of First derivative d2lds2= (d^2l)/(dsigma^2) */ private ArrayRealVector d2lds2() { //d2ldd2 = function(sigma) (4/sigma^4) //-(4/sigma^6)*trigamma((1/sigma^2)), double[] out = new double[size]; for (int i = 0; i < size; i++) { out[i] = (4 / (tempV2.getEntry(i) * tempV2.getEntry(i))) - (4 / (tempV2.getEntry(i) * tempV2.getEntry(i) * tempV2.getEntry(i))) * Gamma.trigamma(1 / (tempV2.getEntry(i))); } return new ArrayRealVector(out, false); }
From source file:gamlss.distributions.TF.java
/** Second derivative d2ldn2= (d^2l)/(dnu^2), * where l - log-likelihood function./*from ww w . j a va2 s . c o m*/ * @param y - vector of values of response variable * @return a vector of second derivative d2ldn2= (d^2l)/(dnu^2) */ private ArrayRealVector d2ldn2(final ArrayRealVector y) { double[] out = new double[size]; for (int i = 0; i < size; i++) { //d2ldv2 <- trigamma(v3)-trigamma(v2)+(2*(nu+5))/(nu*(nu+1)*(nu+3)) out[i] = Gamma.trigamma(v3[i]) - Gamma.trigamma(v2[i]) + (2.0 * (nuV.getEntry(i) + 5)) / (nuV.getEntry(i) * (nuV.getEntry(i) + 1) * (nuV.getEntry(i) + 3)); //d2ldv2 <- d2ldv2/4 out[i] = out[i] / 4.0; //d2ldv2 <- ifelse(d2ldv2 < -1e-15, d2ldv2,-1e-15) if (out[i] > -1e-15) { out[i] = -1e-15; } } return new ArrayRealVector(out, false); }
From source file:gamlss.distributions.TF2.java
/** Second derivative d2ldn2= (d^2l)/(dnu^2), * where l - log-likelihood function.//from w ww . j av a 2s . co m * @param y - vector of values of response variable * @return a vector of second derivative d2ldn2= (d^2l)/(dnu^2) */ private ArrayRealVector d2ldn2(final ArrayRealVector y) { //d2ldv2 <- d2ldv2 + (ds1dv^2)*TF()$d2ldd2(sigma1, nu) tf.setDistributionParameter(DistributionSettings.MU, muV); tf.setDistributionParameter(DistributionSettings.SIGMA, sigma1); tf.setDistributionParameter(DistributionSettings.NU, nuV); tf.firstDerivative(DistributionSettings.SIGMA, y); tempV = tf.secondDerivative(DistributionSettings.SIGMA, y); double[] out = new double[size]; for (int i = 0; i < size; i++) { //v2 <- nu/2 //v3 <-(nu+1)/2 //d2ldv2 <- trigamma(v3)-trigamma(v2)+(2*(nu+5))/(nu*(nu+1)*(nu+3)) //d2ldv2 <- d2ldv2/4 out[i] = (Gamma.trigamma((nuV.getEntry(i) + 1) / 2.0) - Gamma.trigamma(nuV.getEntry(i) / 2.0) + (2.0 * (nuV.getEntry(i) + 5)) / (nuV.getEntry(i) * (nuV.getEntry(i) + 1) * (nuV.getEntry(i) + 3))) / 4.0; out[i] = out[i] + (ds1dv[i] * ds1dv[i]) * tempV.getEntry(i); } sigma1 = null; ds1dv = null; muV = null; sigmaV = null; nuV = null; return new ArrayRealVector(out, false); }
From source file:gamlss.distributions.PE.java
/** Second derivative d2ldn2= (d^2l)/(dnu^2), * where l - log-likelihood function.// w w w . j av a 2 s. c o m * @param y - vector of values of response variable * @return a vector of second derivative d2ldn2= (d^2l)/(dnu^2) */ private ArrayRealVector d2ldn2(final ArrayRealVector y) { double[] out = new double[size]; for (int i = 0; i < size; i++) { //dlogc.dv <- (1/(2*nu^2))*(2*log(2)-digamma(1/nu)+3*digamma(3/nu)) final double dlogcDv = (1 / (2 * nuV.getEntry(i) * nuV.getEntry(i))) * (2 * FastMath.log(2) - Gamma.digamma(1 / nuV.getEntry(i)) + 3 * Gamma.digamma(3 / nuV.getEntry(i))); //p <- (1+nu)/nu final double p = (1 + nuV.getEntry(i)) / nuV.getEntry(i); //part1 <- p*trigamma(p)+2*(digamma(p))^2 final double part1 = p * Gamma.trigamma(p) + 2 * FastMath.pow(Gamma.digamma(p), 2); //part2 <- digamma(p)*(log(2)+3-3*digamma(3/nu)-nu) final double part2 = Gamma.digamma(p) * (FastMath.log(2) + 3 - 3 * Gamma.digamma(3 / nuV.getEntry(i)) - nuV.getEntry(i)); //part3 <- -3*(digamma(3/nu))*(1+log(2)) final double part3 = -3 * (Gamma.digamma(3 / nuV.getEntry(i))) * (1 + FastMath.log(2)); //part4 <- -(nu+log(2))*log(2) final double part4 = -(nuV.getEntry(i) + FastMath.log(2)) * FastMath.log(2); //part5 <- -nu+(nu^4)*(dlogc.dv)^2 final double part5 = -nuV.getEntry(i) + FastMath.pow(nuV.getEntry(i), 4) * dlogcDv * dlogcDv; //d2ldv2 <- part1+part2+part3+part4+part5 out[i] = part1 + part2 + part3 + part4 + part5; //d2ldv2 <- -d2ldv2/nu^3 out[i] = -out[i] / (FastMath.pow(nuV.getEntry(i), 3)); //d2ldv2 <- ifelse(d2ldv2 < -1e-15, d2ldv2,-1e-15) if (out[i] > -1e-15) { out[i] = -1e-15; } } muV = null; sigmaV = null; nuV = null; return new ArrayRealVector(out, false); }