List of usage examples for org.apache.commons.math3.special Beta logBeta
public static double logBeta(final double p, final double q)
From source file:gamlss.distributions.SST.java
/** * Set m1, m2, s1, mu1 arrays./*from w ww. j a v a 2 s.c o m*/ * @param y - response variable */ private void setInterimArrays(final ArrayRealVector y) { muV = distributionParameters.get(DistributionSettings.MU); sigmaV = distributionParameters.get(DistributionSettings.SIGMA); nuV = distributionParameters.get(DistributionSettings.NU); tauV = distributionParameters.get(DistributionSettings.TAU); size = y.getDimension(); m1 = new double[size]; m2 = new double[size]; s1 = new double[size]; mu1 = new double[size]; sigma1 = new double[size]; for (int i = 0; i < size; i++) { //m1 <- (2*tau^(1/2)*(nu^2-1))/((tau-1)*beta(1/2, tau/2)*nu) m1[i] = (2 * (FastMath.sqrt(tauV.getEntry(i))) * (nuV.getEntry(i) * nuV.getEntry(i) - 1)) / ((tauV.getEntry(i) - 1) * FastMath.exp(Beta.logBeta(0.5, tauV.getEntry(i) / 2)) * nuV.getEntry(i)); //m2 <- (tau*(nu^3+(1/nu^3)))/((tau-2)*(nu+(1/nu))) m2[i] = (tauV.getEntry(i) * (nuV.getEntry(i) * nuV.getEntry(i) * nuV.getEntry(i) + (1 / (nuV.getEntry(i) * nuV.getEntry(i) * nuV.getEntry(i))))) / ((tauV.getEntry(i) - 2) * (nuV.getEntry(i) + (1 / nuV.getEntry(i)))); //s1 <- sqrt(m2-m1^2) s1[i] = FastMath.sqrt(m2[i] - m1[i] * m1[i]); //mu1 <- mu- ((sigma*m1)/s1) mu1[i] = muV.getEntry(i) - ((sigmaV.getEntry(i) * m1[i]) / s1[i]); //sigma1 <- sigma/s1 sigma1[i] = sigmaV.getEntry(i) / s1[i]; } }
From source file:gamlss.distributions.SST.java
/** First derivative dldn = dl/dnu, where l - log-likelihood function. * @param y - vector of values of response variable * @return a vector of First derivative dldn = dl/dnu *///from w w w.j a v a2s . co m public final ArrayRealVector dldn(final ArrayRealVector y) { st3.setDistributionParameter(DistributionSettings.MU, new ArrayRealVector(mu1, false)); st3.setDistributionParameter(DistributionSettings.SIGMA, new ArrayRealVector(sigma1, false)); st3.setDistributionParameter(DistributionSettings.NU, nuV); st3.setDistributionParameter(DistributionSettings.TAU, tauV); //dl1dmu1 <- ST3()$dldm(y, mu1, sigma1, nu, tau) //dl1dd1 <- ST3()$dldd(y, mu1, sigma1, nu, tau) //dl1dv <- ST3()$dldv(y, mu1, sigma1, nu, tau) double[] dl1dmu1 = st3.firstDerivative(DistributionSettings.MU, y).getDataRef(); double[] dl1dd1 = st3.firstDerivative(DistributionSettings.SIGMA, y).getDataRef(); double[] dl1dv = st3.firstDerivative(DistributionSettings.NU, y).getDataRef(); dldn = new double[size]; for (int i = 0; i < size; i++) { //dmu1dm1 <- -sigma/s1 final double dmu1dm1 = -(sigmaV.getEntry(i) / s1[i]); //dmu1ds1 <- (sigma*m1)/(s1^2) final double dmu1ds1 = (sigmaV.getEntry(i) * m1[i]) / (s1[i] * s1[i]); //dd1ds1 <- -sigma/(s1^2) final double dd1ds1 = -(sigmaV.getEntry(i) / (s1[i] * s1[i])); //dm1dv <- ((2*tau^(1/2))/((tau-1) //*beta(1/2, tau/2)))*((nu^2+1)/(nu^2)) final double dm1dv = ((2 * FastMath.sqrt(tauV.getEntry(i))) / ((tauV.getEntry(i) - 1) * FastMath.exp(Beta.logBeta(0.5, tauV.getEntry(i) / 2)))) * ((nuV.getEntry(i) * nuV.getEntry(i) + 1) / (nuV.getEntry(i) * nuV.getEntry(i))); //dm2dv <- m2*((6*nu^5/(nu^6+1))-(2/nu)-(2*nu/(nu^2+1))) final double dm2dv = m2[i] * ((6.0 * FastMath.pow(nuV.getEntry(i), 5) / (FastMath.pow(nuV.getEntry(i), 6) + 1)) - (2.0 / nuV.getEntry(i)) - (2.0 * nuV.getEntry(i) / (nuV.getEntry(i) * nuV.getEntry(i) + 1.0))); //ds1dv <- (dm2dv - 2*m1*dm1dv)/(2*s1) final double ds1dv = (dm2dv - 2.0 * m1[i] * dm1dv) / (2 * s1[i]); //dldv <- dl1dmu1*dmu1dm1*dm1dv //+ dl1dmu1*dmu1ds1*ds1dv + dl1dd1*dd1ds1*ds1dv + dl1dv dldn[i] = dl1dmu1[i] * dmu1dm1 * dm1dv + dl1dmu1[i] * dmu1ds1 * ds1dv + dl1dd1[i] * dd1ds1 * ds1dv + dl1dv[i]; } m1 = null; m2 = null; s1 = null; mu1 = null; sigma1 = null; return new ArrayRealVector(dldn, false); }
From source file:gamlss.distributions.SST.java
/** Computes the probability density function (PDF) of this * distribution evaluated at the specified point x. * @param x - value of response variable * @param mu - value of mu distribution parameter * @param sigma - value of sigma distribution parameter * @param nu - value of nu distribution parameter * @param tau - value of tau distribution parameter * @param isLog - logical, whether to take log of the function or not * @return value of probability density function *//*from www . j av a 2 s . c o m*/ public final double dSST(final double x, final double mu, final double sigma, final double nu, final double tau, final boolean isLog) { //m1 <- (2*tau^(1/2)*(nu^2-1))/((tau-1)*beta(1/2, tau/2)*nu) final double m1 = (2 * (FastMath.sqrt(tau)) * (nu * nu - 1)) / ((tau - 1) * FastMath.exp(Beta.logBeta(0.5, tau / 2)) * nu); //m2 <- (tau*(nu^3+(1/nu^3)))/((tau-2)*(nu+(1/nu))) final double m2 = (tau * (nu * nu * nu + (1 / (nu * nu * nu)))) / ((tau - 2) * (nu + (1 / nu))); //s1 <- sqrt(m2-m1^2) final double s1 = FastMath.sqrt(m2 - m1 * m1); //mu1 <- mu- ((sigma*m1)/s1) final double mu1 = mu - ((sigma * m1) / s1); //sigma1 <- sigma/s1 final double sigma1 = sigma / s1; return st3.dST3(x, mu1, sigma1, nu, tau, isLog); }
From source file:gamlss.distributions.SST.java
/** Computes the cumulative distribution * function P(X <= q) for a random variable X . * whose values are distributed according to this distribution * @param q - value of quantile/* ww w.j a v a2 s. co m*/ * @param mu - value of mu distribution parameter * @param sigma - value of sigma distribution parameter * @param nu - value of nu distribution parameter * @param tau - value of tau distribution parameter * @param lowerTail - logical, if TRUE (default), probabilities * are P[X <= x] otherwise, P[X > x]. * @param isLog - logical, if TRUE, probabilities p are given as log(p) * @return value of cumulative probability function values P(X <= q) */ public final double pSST(final double q, final double mu, final double sigma, final double nu, final double tau, final boolean lowerTail, final boolean isLog) { //m1 <- (2*tau^(1/2)*(nu^2-1))/((tau-1)*beta(1/2, tau/2)*nu) final double m1 = (2 * (FastMath.sqrt(tau)) * (nu * nu - 1)) / ((tau - 1) * FastMath.exp(Beta.logBeta(0.5, tau / 2)) * nu); //m2 <- (tau*(nu^3+(1/nu^3)))/((tau-2)*(nu+(1/nu))) final double m2 = (tau * (nu * nu * nu + (1 / (nu * nu * nu)))) / ((tau - 2) * (nu + (1 / nu))); //s1 <- sqrt(m2-m1^2) final double s1 = FastMath.sqrt(m2 - m1 * m1); //mu1 <- mu- ((sigma*m1)/s1) final double mu1 = mu - ((sigma * m1) / s1); //sigma1 <- sigma/s1 final double sigma1 = sigma / s1; return st3.pST3(q, mu1, sigma1, nu, tau, lowerTail, isLog); }
From source file:gamlss.distributions.SST.java
/** Computes the quantile (inverse cumulative probability) * function of this distribution./* www. java2 s . c om*/ * @param p - value of cumulative probability * @param mu - value of mu distribution parameters * @param sigma - value of sigma distribution parameters * @param nu - value of nu distribution parameters * @param tau - value of tau distribution parameters * @param lowerTail - logical; if TRUE (default), probabilities * are P[X <= x] otherwise, P[X > x] * @param isLog - logical; if TRUE, probabilities p are given as log(p). * @return value of quantile function */ public final double qSST(final double p, final double mu, final double sigma, final double nu, final double tau, final boolean lowerTail, final boolean isLog) { //if (any(tau <= 2)) stop(paste("tau must be greater than 2")) if (tau <= 2) { System.err.println("tau must be greater than 2"); return -1.0; } //m1 <- (2*tau^(1/2)*(nu^2-1))/((tau-1)*beta(1/2, tau/2)*nu) final double m1 = (2 * (FastMath.sqrt(tau)) * (nu * nu - 1)) / ((tau - 1) * FastMath.exp(Beta.logBeta(0.5, tau / 2)) * nu); //m2 <- (tau*(nu^3+(1/nu^3)))/((tau-2)*(nu+(1/nu))) final double m2 = (tau * (nu * nu * nu + (1 / (nu * nu * nu)))) / ((tau - 2) * (nu + (1 / nu))); //s1 <- sqrt(m2-m1^2) final double s1 = FastMath.sqrt(m2 - m1 * m1); //mu1 <- mu- ((sigma*m1)/s1) final double mu1 = mu - ((sigma * m1) / s1); //sigma1 <- sigma/s1 final double sigma1 = sigma / s1; return st3.qST3(p, mu1, sigma1, nu, tau, lowerTail, isLog); }
From source file:nl.systemsgenetics.cellTypeSpecificAlleleSpecificExpression.LikelihoodFunctions.java
static double BetaBinomLogLik(double sigma, double alpha, double beta, int asRef, int asAlt) { int AS1 = asRef; int AS2 = asAlt; double hetp = GlobalVariables.hetProb; double error = GlobalVariables.seqError; double a;//from w ww .j ava 2 s . co m double b; a = Math.exp((Math.log(alpha) - Math.log(alpha + beta)) + Math.log((1.0 / Math.pow(sigma, 2)) - 1.0)); b = Math.exp((Math.log(beta) - Math.log(alpha + beta)) + Math.log((1.0 / Math.pow(sigma, 2)) - 1.0)); double part1 = 0.0; part1 += Beta.logBeta(AS1 + a, AS2 + b); part1 -= Beta.logBeta(a, b); if (GlobalVariables.hetProb == 1.0) { return -1.0 * part1; } double e1 = Math.log(error) * AS1 + Math.log(1.0 - error) * AS2; double e2 = Math.log(error) * AS2 + Math.log(1.0 - error) * AS1; if (GlobalVariables.hetProb == 0.0) { return -1.0 * addlogs(e1, e2); } return -1.0 * addlogs(Math.log(hetp) + part1, Math.log(1 - hetp) + addlogs(e1, e2)); }
From source file:uk.ac.diamond.scisoft.analysis.fitting.functions.PearsonVII.java
@Override protected void calcCachedParameters() { pos = getParameterValue(POSN);/*from w w w . j av a 2 s . co m*/ power = getParameterValue(POWER); halfwp = 0.5 * getParameterValue(FWHM) / Math.sqrt(Math.pow(2, 1. / power) - 1); double beta = Math.exp(Beta.logBeta(power - 0.5, 0.5)); height = getParameterValue(AREA) / (beta * halfwp); setDirty(false); }