List of usage examples for org.apache.commons.math3.special Gamma digamma
public static double digamma(double x)
Computes the digamma function of x.
This is an independently written implementation of the algorithm described in Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976.
Some of the constants have been changed to increase accuracy at the moderate expense of run-time.
From source file:gamlss.distributions.BCPE.java
/** First derivative dldtau = dl/dtau, * where l - log-likelihood function.//from w ww.jav a2 s. com * @param y - vector of values of response variable * @return a vector of First derivative dldtau = dl/dtau */ public final ArrayRealVector dldt(final ArrayRealVector y) { dldt = new double[size]; for (int i = 0; i < size; i++) { //j <- (log(F.T(1/(sigma*abs(nu)),tau+0.001)) //-log(F.T(1/(sigma*abs(nu)),tau)))/0.001 final double j = (FastMath .log(f2T(1 / (sigmaV.getEntry(i) * FastMath.abs(nuV.getEntry(i))), tauV.getEntry(i) + 0.001, i)) - FastMath.log( f2T(1 / (sigmaV.getEntry(i) * FastMath.abs(nuV.getEntry(i))), tauV.getEntry(i), i))) / 0.001; //dldt <- (1/tau)-0.5*(log(abs(z/c)))* //(abs(z/c))^tau+(1/tau^2)*(log(2)+digamma(1/tau))+((tau/2)* //((abs(z/c))^tau)-1)*dlogc.dt-j dldt[i] = (1 / tauV.getEntry(i)) - 0.5 * (FastMath.log(FastMath.abs(z[i] / c[i]))) * (FastMath.pow(FastMath.abs(z[i] / c[i]), tauV.getEntry(i))) + (1 / (tauV.getEntry(i) * tauV.getEntry(i))) * (FastMath.log(2) + Gamma.digamma(1 / tauV.getEntry(i))) + ((tauV.getEntry(i) / 2) * (FastMath.pow(FastMath.abs(z[i] / c[i]), tauV.getEntry(i))) - 1) * dlogcDt[i] - j; } c = null; logC = null; z = null; return new ArrayRealVector(dldt, false); }
From source file:gamlss.distributions.PE.java
/** Second cross derivative of likelihood function * in respect to sigma and nu (d2ldmdd = d2l/dsigma*dnu). * @param y - vector of values of response variable * @return a vector of Second cross derivative *///from www . j av a2 s .com private ArrayRealVector d2ldsdn(final ArrayRealVector y) { ArrayRealVector sigmaT = distributionParameters.get(DistributionSettings.SIGMA); ArrayRealVector nuT = distributionParameters.get(DistributionSettings.NU); double[] out = new double[size]; for (int i = 0; i < size; i++) { //d2ldddv = function(y,mu,sigma,nu) (1/(2*sigma)) //*((3/nu)*(digamma(1/nu)-digamma(3/nu))+2+(2/nu)) out[i] = (1 / (2 * sigmaT.getEntry(i))) * ((3 / nuT.getEntry(i)) * (Gamma.digamma(1 / nuT.getEntry(i)) - Gamma.digamma(3 / nuT.getEntry(i))) + 2 + (2 / nuT.getEntry(i))); } return new ArrayRealVector(out, false); }
From source file:gamlss.distributions.BCPE.java
/** Second derivative d2ldt2= (d^2l)/(dtau^2), * where l - log-likelihood function.//from ww w . ja va2 s. co m * @param y - vector of values of response variable * @return a vector of second derivative d2ldt2= (d^2l)/(dtau^2) */ private ArrayRealVector d2ldt2(final ArrayRealVector y) { double[] out = new double[size]; for (int i = 0; i < size; i++) { //p <- (tau+1)/tau final double p = (tauV.getEntry(i) + 1) / tauV.getEntry(i); //part1 <- p*trigamma(p)+2*(digamma(p))^2 final double part1 = p * Gamma.trigamma(p) + 2 * (Gamma.digamma(p)) * (Gamma.digamma(p)); //part2 <- digamma(p)*(log(2)+3-3*digamma(3/tau)-tau) final double part2 = Gamma.digamma(p) * (FastMath.log(2) + 3 - 3 * Gamma.digamma(3 / tauV.getEntry(i)) - tauV.getEntry(i)); //part3 <- -3*(digamma(3/tau))*(1+log(2)) final double part3 = -3 * (Gamma.digamma(3 / tauV.getEntry(i))) * (1 + FastMath.log(2)); //part4 <- -(tau+log(2))*log(2) final double part4 = -(tauV.getEntry(i) + FastMath.log(2)) * FastMath.log(2); //part5 <- -tau+(tau^4)*(dlogc.dt)^2 final double part5 = -tauV.getEntry(i) + (FastMath.pow(tauV.getEntry(i), 4)) * (dlogcDt[i]) * (dlogcDt[i]); //d2ldt2 <- part1+part2+part3+part4+part5 out[i] = part1 + part2 + part3 + part4 + part5; //d2ldt2 <- -d2ldt2/tau^3 out[i] = -out[i] / (tauV.getEntry(i) * tauV.getEntry(i) * tauV.getEntry(i)); //d2ldt2 <- ifelse(d2ldt2 < -1e-15, d2ldt2,-1e-15) if (out[i] > -1e-15) { out[i] = -1e-15; } } muV = null; sigmaV = null; nuV = null; tauV = null; return new ArrayRealVector(out, false); }
From source file:gamlss.distributions.BCPE.java
/** Second cross derivative of likelihood function * in respect to mu and tau (d2ldmdd = d2l/dmu*dtau). * @param y - vector of values of response variable * @return a vector of Second cross derivative *///from w w w .j a va2 s .c o m private ArrayRealVector d2ldmdt(final ArrayRealVector y) { ArrayRealVector mu = distributionParameters.get(DistributionSettings.MU); ArrayRealVector sigma = distributionParameters.get(DistributionSettings.SIGMA); ArrayRealVector nu = distributionParameters.get(DistributionSettings.NU); ArrayRealVector tau = distributionParameters.get(DistributionSettings.TAU); size = y.getDimension(); double[] out = new double[size]; for (int i = 0; i < size; i++) { //d2ldmdt = (nu/(mu*tau))*(1+tau+(3/2) //*(digamma(1/tau)-digamma(3/tau))) out[i] = (nu.getEntry(i) / (mu.getEntry(i) * tau.getEntry(i))) * (1 + tau.getEntry(i) + (3 / 2) * (Gamma.digamma(1 / tau.getEntry(i)) - Gamma.digamma(3 / tau.getEntry(i)))); } return new ArrayRealVector(out, false); }
From source file:gamlss.distributions.BCPE.java
/** Second cross derivative of likelihood function * in respect to sigma and tau (d2ldmdd = d2l/dsigma*dtau). * @param y - vector of values of response variable * @return a vector of Second cross derivative *///w ww .j a v a 2 s . co m private ArrayRealVector d2ldsdt(final ArrayRealVector y) { ArrayRealVector mu = distributionParameters.get(DistributionSettings.MU); ArrayRealVector sigma = distributionParameters.get(DistributionSettings.SIGMA); ArrayRealVector nu = distributionParameters.get(DistributionSettings.NU); ArrayRealVector tau = distributionParameters.get(DistributionSettings.TAU); size = y.getDimension(); double[] out = new double[size]; for (int i = 0; i < size; i++) { //d2ldddt = (1/(sigma*tau))*(1+tau+(3/2) //*(digamma(1/tau)-digamma(3/tau))) out[i] = (1 / (sigma.getEntry(i) * tau.getEntry(i))) * (1 + tau.getEntry(i) + (3 / 2) * (Gamma.digamma(1 / tau.getEntry(i)) - Gamma.digamma(3 / tau.getEntry(i)))); } return new ArrayRealVector(out, false); }
From source file:gamlss.distributions.BCPE.java
/** Second cross derivative of likelihood function * in respect to nu and tau (d2ldmdd = d2l/dnu*dtau). * @param y - vector of values of response variable * @return a vector of Second cross derivative *//*from w w w . jav a 2 s . c o m*/ private ArrayRealVector d2ldndt(final ArrayRealVector y) { ArrayRealVector mu = distributionParameters.get(DistributionSettings.MU); ArrayRealVector sigma = distributionParameters.get(DistributionSettings.SIGMA); ArrayRealVector nu = distributionParameters.get(DistributionSettings.NU); ArrayRealVector tau = distributionParameters.get(DistributionSettings.TAU); size = y.getDimension(); double[] out = new double[size]; for (int i = 0; i < size; i++) { //d2ldvdt <- (((sigma^2)*nu)/(2*tau)) //*(1+(tau/3)+0.5*(digamma(1/tau)-digamma(3/tau))) out[i] = (((sigma.getEntry(i) * sigma.getEntry(i)) * nu.getEntry(i)) / (2 * tau.getEntry(i))) * (1 + (tau.getEntry(i) / 3) + 0.5 * (Gamma.digamma(1 / tau.getEntry(i)) - Gamma.digamma(3 / tau.getEntry(i)))); } return new ArrayRealVector(out, false); }
From source file:opennlp.tools.fca.Measures.java
public void probability() { //the probability of B being closed for (int i = 0; i < cl.conceptList.size(); i++) { ArrayList<Integer> intent = cl.conceptList.get(i).getIntent(); // out of concept intent double pB = intentProbability(intent); ArrayList<Integer> outOfIntent = new ArrayList<Integer>(); ArrayList<Double> outOfIntentAttrProb = new ArrayList<Double>(); for (int j = 0; j < cl.attributeCount; j++) { outOfIntent.add(j);/*from ww w.java 2s. c o m*/ } for (int j = intent.size() - 1; j >= 0; j--) { outOfIntent.remove(intent.get(j)); } for (int j = 0; j < outOfIntent.size(); j++) { outOfIntentAttrProb.add(attributeProbability(outOfIntent.get(j))); } double prob = 0, mult = 1, mult1 = 1; int n = cl.objectCount; for (int k = 0; k <= n; k++) { mult = 1; mult1 = 1; for (int j = 0; j < outOfIntentAttrProb.size(); j++) { mult *= (1 - Math.pow(outOfIntentAttrProb.get(j), k)); } mult1 = Math.pow(pB, k) * Math.pow(1 - pB, n - k); prob += mult1 * mult * Gamma.digamma(n + 1) / Gamma.digamma(k + 1) / Gamma.digamma(n - k + 1); } cl.conceptList.get(i).probability = prob; } }
From source file:org.broadinstitute.gatk.utils.Dirichlet.java
public double[] effectiveMultinomialWeights() { final double digammaOfSum = Gamma.digamma(MathUtils.sum(alpha)); return MathUtils.applyToArray(alpha, a -> Math.exp(Gamma.digamma(a) - digammaOfSum)); }
From source file:org.broadinstitute.gatk.utils.Dirichlet.java
public double[] effectiveLog10MultinomialWeights() { final double digammaOfSum = Gamma.digamma(MathUtils.sum(alpha)); return MathUtils.applyToArray(alpha, a -> (Gamma.digamma(a) - digammaOfSum) * MathUtils.LOG10_OF_E); }
From source file:org.lenskit.pf.HPFModelProvider.java
@Override public HPFModel get() { final int userNum = ratings.getUserIndex().size(); final int itemNum = ratings.getItemIndex().size(); final int featureCount = hyperParameters.getFeatureCount(); final double a = hyperParameters.getUserWeightShpPrior(); final double aPrime = hyperParameters.getUserActivityShpPrior(); final double bPrime = hyperParameters.getUserActivityPriorMean(); final double c = hyperParameters.getItemWeightShpPrior(); final double cPrime = hyperParameters.getItemActivityShpPrior(); final double dPrime = hyperParameters.getItemActivityPriorMean(); final double kappaShpU = aPrime + featureCount * a; final double tauShpI = cPrime + featureCount * c; RealMatrix gammaShp = MatrixUtils.createRealMatrix(userNum, featureCount); RealMatrix gammaRte = MatrixUtils.createRealMatrix(userNum, featureCount); RealVector kappaShp = MatrixUtils.createRealVector(new double[userNum]); RealVector kappaRte = MatrixUtils.createRealVector(new double[userNum]); RealMatrix lambdaShp = MatrixUtils.createRealMatrix(itemNum, featureCount); RealMatrix lambdaRte = MatrixUtils.createRealMatrix(itemNum, featureCount); RealVector tauShp = MatrixUtils.createRealVector(new double[itemNum]); RealVector tauRte = MatrixUtils.createRealVector(new double[itemNum]); RealMatrix gammaShpNext = MatrixUtils.createRealMatrix(userNum, featureCount); RealMatrix lambdaShpNext = MatrixUtils.createRealMatrix(itemNum, featureCount); gammaShpNext = gammaShpNext.scalarAdd(a); lambdaShpNext = lambdaShpNext.scalarAdd(c); RealVector phiUI = MatrixUtils.createRealVector(new double[featureCount]); initialize(gammaShp, gammaRte, kappaRte, kappaShp, lambdaShp, lambdaRte, tauRte, tauShp); logger.info("initialization finished"); final List<RatingMatrixEntry> train = ratings.getTrainRatings(); final List<RatingMatrixEntry> validation = ratings.getValidationRatings(); double avgPLLPre = Double.MAX_VALUE; double avgPLLCurr = 0.0; double diffPLL = 1.0; int iterCount = 1; while (iterCount < maxIterCount && diffPLL > threshold) { // update phi Iterator<RatingMatrixEntry> allUIPairs = train.iterator(); while (allUIPairs.hasNext()) { RatingMatrixEntry entry = allUIPairs.next(); int item = entry.getItemIndex(); int user = entry.getUserIndex(); double ratingUI = entry.getValue(); if (ratingUI <= 0) { continue; }/*from ww w . j a v a2 s. co m*/ for (int k = 0; k < featureCount; k++) { double gammaShpUK = gammaShp.getEntry(user, k); double gammaRteUK = gammaRte.getEntry(user, k); double lambdaShpIK = lambdaShp.getEntry(item, k); double lambdaRteIK = lambdaRte.getEntry(item, k); double phiUIK = Gamma.digamma(gammaShpUK) - Math.log(gammaRteUK) + Gamma.digamma(lambdaShpIK) - Math.log(lambdaRteIK); phiUI.setEntry(k, phiUIK); } logNormalize(phiUI); if (ratingUI > 1) { phiUI.mapMultiplyToSelf(ratingUI); } for (int k = 0; k < featureCount; k++) { double value = phiUI.getEntry(k); gammaShpNext.addToEntry(user, k, value); lambdaShpNext.addToEntry(item, k, value); } } logger.info("iteration {} first phrase update finished", iterCount); RealVector gammaRteSecondTerm = MatrixUtils.createRealVector(new double[featureCount]); for (int k = 0; k < featureCount; k++) { double gammaRteUK = 0.0; for (int item = 0; item < itemNum; item++) { gammaRteUK += lambdaShp.getEntry(item, k) / lambdaRte.getEntry(item, k); } gammaRteSecondTerm.setEntry(k, gammaRteUK); } // update user parameters double kappaRteFirstTerm = aPrime / bPrime; for (int user = 0; user < userNum; user++) { double gammaRteUKFirstTerm = kappaShp.getEntry(user) / kappaRte.getEntry(user); double kappaRteU = 0.0; for (int k = 0; k < featureCount; k++) { double gammaShpUK = gammaShpNext.getEntry(user, k); gammaShp.setEntry(user, k, gammaShpUK); gammaShpNext.setEntry(user, k, a); double gammaRteUK = gammaRteSecondTerm.getEntry(k); gammaRteUK += gammaRteUKFirstTerm; gammaRte.setEntry(user, k, gammaRteUK); kappaRteU += gammaShpUK / gammaRteUK; } kappaRteU += kappaRteFirstTerm; kappaRte.setEntry(user, kappaRteU); } logger.info("iteration {} second phrase update finished", iterCount); RealVector lambdaRteSecondTerm = MatrixUtils.createRealVector(new double[featureCount]); for (int k = 0; k < featureCount; k++) { double lambdaRteIK = 0.0; for (int user = 0; user < userNum; user++) { lambdaRteIK += gammaShp.getEntry(user, k) / gammaRte.getEntry(user, k); } lambdaRteSecondTerm.setEntry(k, lambdaRteIK); } // update item parameters double tauRteFirstTerm = cPrime / dPrime; for (int item = 0; item < itemNum; item++) { double lambdaRteFirstTerm = tauShp.getEntry(item) / tauRte.getEntry(item); double tauRteI = 0.0; for (int k = 0; k < featureCount; k++) { double lambdaShpIK = lambdaShpNext.getEntry(item, k); lambdaShp.setEntry(item, k, lambdaShpIK); lambdaShpNext.setEntry(item, k, c); double lambdaRteIK = lambdaRteSecondTerm.getEntry(k); // plus first term lambdaRteIK += lambdaRteFirstTerm; lambdaRte.setEntry(item, k, lambdaRteIK); // update tauRteI second term tauRteI += lambdaShpIK / lambdaRteIK; } tauRteI += tauRteFirstTerm; tauRte.setEntry(item, tauRteI); } logger.info("iteration {} third phrase update finished", iterCount); // compute average predictive log likelihood of validation data per {@code iterationfrequency} iterations if (iterCount == 1) { for (int user = 0; user < userNum; user++) { kappaShp.setEntry(user, kappaShpU); } for (int item = 0; item < itemNum; item++) { tauShp.setEntry(item, tauShpI); } } if ((iterCount % iterationFrequency) == 0) { Iterator<RatingMatrixEntry> valIter = validation.iterator(); avgPLLCurr = 0.0; while (valIter.hasNext()) { RatingMatrixEntry ratingEntry = valIter.next(); int user = ratingEntry.getUserIndex(); int item = ratingEntry.getItemIndex(); double rating = ratingEntry.getValue(); double eThetaBeta = 0.0; for (int k = 0; k < featureCount; k++) { double eThetaUK = gammaShp.getEntry(user, k) / gammaRte.getEntry(user, k); double eBetaIK = lambdaShp.getEntry(item, k) / lambdaRte.getEntry(item, k); eThetaBeta += eThetaUK * eBetaIK; } double pLL = 0.0; if (isProbPredition) { pLL = (rating == 0) ? (-eThetaBeta) : Math.log(1 - Math.exp(-eThetaBeta)); } else { pLL = rating * Math.log(eThetaBeta) - eThetaBeta - Gamma.logGamma(rating + 1); } avgPLLCurr += pLL; } avgPLLCurr = avgPLLCurr / validation.size(); diffPLL = Math.abs((avgPLLCurr - avgPLLPre) / avgPLLPre); avgPLLPre = avgPLLCurr; logger.info("iteration {} with current average predictive log likelihood {} and the change is {}", iterCount, avgPLLCurr, diffPLL); } iterCount++; } // construct feature matrix used by HPFModel RealMatrix eTheta = MatrixUtils.createRealMatrix(userNum, featureCount); RealMatrix eBeta = MatrixUtils.createRealMatrix(itemNum, featureCount); for (int user = 0; user < userNum; user++) { RealVector gammaShpU = gammaShp.getRowVector(user); RealVector gammaRteU = gammaRte.getRowVector(user); RealVector eThetaU = gammaShpU.ebeDivide(gammaRteU); eTheta.setRowVector(user, eThetaU); logger.info("Training user {} features finished", user); } for (int item = 0; item < itemNum; item++) { RealVector lambdaShpI = lambdaShp.getRowVector(item); RealVector lambdaRteI = lambdaRte.getRowVector(item); RealVector eBetaI = lambdaShpI.ebeDivide(lambdaRteI); eBeta.setRowVector(item, eBetaI); logger.info("Training item {} features finished", item); } KeyIndex uidx = ratings.getUserIndex(); KeyIndex iidx = ratings.getItemIndex(); return new HPFModel(eTheta, eBeta, uidx, iidx); }