List of usage examples for org.apache.commons.math.special Gamma logGamma
public static double logGamma(double x)
From source file:net.malariagen.gatk.math.SaddlePointExpansion.java
/** * Compute the error of Stirling's series at the given value. * <p>// w w w .ja v a2 s . co m * References: * <ol> * <li>Eric W. Weisstein. "Stirling's Series." From MathWorld--A Wolfram Web * Resource. <a target="_blank" * href="http://mathworld.wolfram.com/StirlingsSeries.html"> * http://mathworld.wolfram.com/StirlingsSeries.html</a></li> * </ol> * </p> * * @param z the value. * @return the Striling's series error. */ static double getStirlingError(double z) { double ret; if (z < 15.0) { double z2 = 2.0 * z; if (FastMath.floor(z2) == z2) { ret = EXACT_STIRLING_ERRORS[(int) z2]; } else { ret = Gamma.logGamma(z + 1.0) - (z + 0.5) * FastMath.log(z) + z - HALF_LOG_2_PI; } } else { double z2 = z * z; ret = (0.083333333333333333333 - (0.00277777777777777777778 - (0.00079365079365079365079365 - (0.000595238095238095238095238 - 0.0008417508417508417508417508 / z2) / z2) / z2) / z2) / z; } return ret; }
From source file:com.opengamma.analytics.math.statistics.distribution.NonCentralChiSquaredDistribution.java
/** * {@inheritDoc}// w ww . j a v a 2 s .c o m */ @Override public double getCDF(final Double x) { Validate.notNull(x, "x"); if (x < 0) { return 0.0; } if ((_dofOverTwo + _lambdaOverTwo) > 1000) { return getFraserApproxCDF(x); } double regGammaStart = 0; final double halfX = x / 2.0; final double logX = Math.log(halfX); try { regGammaStart = Gamma.regularizedGammaP(_dofOverTwo + _k, halfX); } catch (final org.apache.commons.math.MathException ex) { throw new MathException(ex); } double sum = _pStart * regGammaStart; double oldSum = Double.NEGATIVE_INFINITY; double p = _pStart; double regGamma = regGammaStart; double temp; int i = _k; // first add terms below _k while (i > 0 && Math.abs(sum - oldSum) / sum > _eps) { i--; p *= (i + 1) / _lambdaOverTwo; temp = (_dofOverTwo + i) * logX - halfX - Gamma.logGamma(_dofOverTwo + i + 1); regGamma += Math.exp(temp); oldSum = sum; sum += p * regGamma; } p = _pStart; regGamma = regGammaStart; oldSum = Double.NEGATIVE_INFINITY; i = _k; while (Math.abs(sum - oldSum) / sum > _eps) { i++; p *= _lambdaOverTwo / i; temp = (_dofOverTwo + i - 1) * logX - halfX - Gamma.logGamma(_dofOverTwo + i); regGamma -= Math.exp(temp); oldSum = sum; sum += p * regGamma; } return sum; }
From source file:etymology.util.EtyMath.java
public static double lnGamma(double value) { if (value < 0) { throw new RuntimeException("Value for lnGamma must be > 0."); }//w ww. ja v a 2s . com if (value == 0) { return 0; } if (value == 1.0) { return 0; } if (value > 100000) { return avgLnGamma(value); } int arrayIndex = (int) Math.ceil((value * GAMMA_CACHE_DIVIDER)); if (logGammaValues.size() > arrayIndex) { return logGammaValues.get(arrayIndex); } int size = (logGammaValues.isEmpty()) ? 1 : logGammaValues.size(); double val = 0.0; for (int i = size; i <= arrayIndex; i++) { val = Gamma.logGamma((i / GAMMA_CACHE_DIVIDER)); logGammaValues.add(val); } return val; }
From source file:geogebra.util.MyMath.java
/** * Factorial function of x. If x is an integer value x! is returned, * otherwise gamma(x + 1) will be returned. For x < 0 Double.NaN is * returned./*from www . j a va2s . c o m*/ * @param x * @return factorial */ final public static double factorial(double x) { if (x < 0) return Double.NaN; // bugfix Michael Borcherds 2008-05-04 // big x or floating point x is computed using gamma function if (x < 0 || x > 32 || x - Math.floor(x) > 1E-10) // exp of log(gamma(x+1)) return Math.exp(Gamma.logGamma(x + 1.0)); int n = (int) x; int j; while (factorialTop < n) { j = factorialTop++; factorialTable[factorialTop] = factorialTable[j] * factorialTop; } return factorialTable[n]; }
From source file:geogebra.kernel.AlgoBinomial.java
private double BinomLog(double n, double r) { // exact for n<=37 // also if r<2.8+Math.exp((250-n)/100) && n<59000 // eg Binom2(38,19) is wrong return Math.floor( 0.5 + Math.exp(Gamma.logGamma(n + 1d) - Gamma.logGamma(r + 1) - Gamma.logGamma((n - r) + 1))); }
From source file:com.hmsinc.epicenter.spatial.analysis.SpatialScanGrid.java
/** Analyze the grid to calculate values for C_all and B_all */ public void analyzeGrid() { B_all = 0.0;// w ww .j a v a 2s . co m C_all = 0.0; double sum_q = 0.0; // Iterate through the grid, summing all baseline and count values. We also // calculate the sum of C/B ratios, for possible future use (these would be // necessary to use population values as the baseline) for (int i = 0; i < rows; i++) { for (int j = 0; j < columns; j++) { B_all += grid[i][j].getBaseline(); C_all += grid[i][j].getCount(); // increment sum_q if baseline != 0 (to avoid NaN) if (grid[i][j].getBaseline() > 0.0) { sum_q += grid[i][j].getCount() / grid[i][j].getBaseline(); } numCells++; } } // Disabled because we use expected value as the baseline...in this case, q0 // is hard-coded to 1.0 // q0 = sum_q/numCells; double alpha_all = q0 * B_all; double beta_all = B_all; // calculate log likelihood of the null hypothesis /** Note: Individual terms in this equation are too large, and return * infinite results. So, we reformulate the equation by taking its * natural log and finding its exponential at the end of calculation. * * Thus the original equation (Neill et al, page 3) * * P(D|H0) ~ (beta_all)^alpha_all * Gamma(alpha_all + C_all) * -------------------------------------------------------- * (beta_all + B_all)^(alpha_all + C_all) * Gamma(alpha_all) * * becomes * * log(P(D|H0)) ~ [ ( alpha_all * log(beta_all) ) + logGamma(alpha_all + C_all) ] - * [ ( (alpha_all + C_all) * log(beta_all + B_all) ) + logGamma(alpha_all) ] * */ logP_D_H0 = ((alpha_all * Math.log(beta_all)) + Gamma.logGamma((alpha_all + C_all))) - (((alpha_all + C_all) * (Math.log(beta_all + B_all))) + Gamma.logGamma(alpha_all)); }
From source file:com.hmsinc.epicenter.spatial.analysis.BayesianSpatialScanStatistic.java
/** Calculate the log Likelihood of an outbreak in a given region */ private double calculateLogLikelihood(SpatialScanGrid grid, double mInitial, double mFinal, double mIncrement, double C_in, double C_out, double B_in, double B_out) { int mCounter = 0; double logLikelihoodSum = 0.0; logger.trace("C_in: {}, C_out: {}, B_in: {}, B_out: {}, Q0: {}", new Object[] { C_in, C_out, B_in, B_out, grid.getQ0() }); /** /*w ww .ja v a 2s. co m*/ * This loop exists to account for the unknown value of parameter m. Final * scores are calculated by averaging over the distribution of m from * [mInitial:mFinal] at mIncrement intervals. * */ for (double m = mInitial; m <= mFinal; m += mIncrement) { double alpha_in = m * grid.getQ0() * B_in; double beta_in = B_in; double alpha_out = grid.getQ0() * B_out; double beta_out = B_out; /** We have all of the values we need, now to calculate probabilities */ /** Note: Individual terms in this equation are too large, and return * infinite results. So, we reformulate the equation by taking its * natural log. Eventually probability will be converted back to a * non-log form. * * E.g. log(beta_in^alpha_in) = alpha_in * log(beta_in) * * Thus the original equation * * (beta_in)^alpha_in * Gamma(alpha_in + C_in) (beta_out)^alpha_out * Gamma(alpha_out + C_out) * P(D|H1(S)) ~ ---------------------------------------------------- X --------------------------------------------------------- * (beta_in + B_in)^(alpha_in + C_in) * Gamma(alpha_in) (beta_out + B_out)^(alpha_out + C_out) * Gamma(alpha_out) * * becomes * * log(P(D|H1(S))) ~ [ ( alpha_in * log(beta_in) ) + logGamma(alpha_in + C_in) ] - * [ ( (alpha_in + C_in) * log(beta_in + B_in) ) + logGamma(alpha_inl) ] + * [ ( alpha_out * log(beta_out) ) + logGamma(alpha_out + C_out) ] - * [ ( (alpha_out + C_out) * log(beta_out + B_out) ) + logGamma(alpha_out) ] * */ double logLikelihood = (((alpha_in * Math.log(beta_in)) + Gamma.logGamma((alpha_in + (double) C_in))) - (((alpha_in + (double) C_in) * (Math.log(beta_in + (double) B_in))) + Gamma.logGamma(alpha_in))) + (((alpha_out * Math.log(beta_out)) + Gamma.logGamma((alpha_out + (double) C_out))) - (((alpha_out + (double) C_out) * (Math.log(beta_out + (double) B_out))) + Gamma.logGamma(alpha_out))); logLikelihoodSum += logLikelihood; mCounter++; } //end for loop over m Validate.isTrue(mCounter > 0, "mCounter was 0!"); double logLikelihood = (logLikelihoodSum / ((double) mCounter)); return logLikelihood; }
From source file:org.apache.mahout.clustering.lda.LDAInference.java
/** * diGamma(x) = gamma'(x)/gamma(x)/*www .j a v a2 s . c o m*/ * logGamma(x) = log(gamma(x)) * * ll = log(gamma(smooth*numTop) / smooth^numTop) + * sum_{i < numTop} (smooth - g[i])*(digamma(g[i]) - digamma(|g|)) + log(gamma(g[i]) * Computes the log likelihood of the wordCounts vector, given \phi, \gamma, and \digamma(gamma) * @param wordCounts * @param map * @param phi * @param gamma * @param digammaGamma * @return */ private double computeLikelihood(Vector wordCounts, int[] map, Matrix phi, Vector gamma, Vector digammaGamma) { double ll = 0.0; // log normalizer for q(gamma); ll += Gamma.logGamma(state.getTopicSmoothing() * state.getNumTopics()); ll -= state.getNumTopics() * Gamma.logGamma(state.getTopicSmoothing()); // isNotNaNAssertion(ll); // now for the the rest of q(gamma); for (int k = 0; k < state.getNumTopics(); ++k) { double gammaK = gamma.get(k); ll += (state.getTopicSmoothing() - gammaK) * digammaGamma.getQuick(k); ll += Gamma.logGamma(gammaK); } ll -= Gamma.logGamma(gamma.zSum()); // isNotNaNAssertion(ll); // for each word for (Iterator<Vector.Element> iter = wordCounts.iterateNonZero(); iter.hasNext();) { Vector.Element e = iter.next(); int w = e.index(); double n = e.get(); int mapping = map[w]; // now for each topic: for (int k = 0; k < state.getNumTopics(); k++) { double llPart = 0.0; double phiKMapping = phi.getQuick(k, mapping); llPart += Math.exp(phiKMapping) * (digammaGamma.getQuick(k) - phiKMapping + state.logProbWordGivenTopic(w, k)); ll += llPart * n; // likelihoodAssertion(w, k, llPart); } } // isLessThanOrEqualsZero(ll); return ll; }
From source file:org.earthtime.UPb_Redux.dateInterpretation.WeightedMeanGraphPanel.java
/** * * @param g2d/*from w ww. j a v a2s . co m*/ */ public void paint(Graphics2D g2d) { // setup painting parameters String fractionSortOrder = "name"; //random, weight, date if (getWeightedMeanOptions().containsKey("fractionSortOrder")) { fractionSortOrder = getWeightedMeanOptions().get("fractionSortOrder"); } double rangeX = (getMaxX_Display() - getMinX_Display()); double rangeY = (getMaxY_Display() - getMinY_Display()); g2d.setClip(getLeftMargin(), getTopMargin(), getGraphWidth(), getGraphHeight()); RenderingHints rh = g2d.getRenderingHints(); rh.put(RenderingHints.KEY_ANTIALIASING, RenderingHints.VALUE_ANTIALIAS_ON); rh.put(RenderingHints.KEY_RENDERING, RenderingHints.VALUE_RENDER_QUALITY); g2d.setRenderingHints(rh); // walk the sampleDateInterpretations and produce graphs g2d.setPaint(Color.BLACK); g2d.setStroke(new BasicStroke(2.0f)); g2d.setFont(new Font("SansSerif", Font.BOLD, 10)); double barWidth = 15.0; double barGap = 10.0; double startSamX = 10.0; double saveStartSamX = 0.0; double samSpace = 3.0; for (int i = 0; i < selectedSampleDateModels.length; i++) { for (int j = 1; j < 9; j++) { if (selectedSampleDateModels[i][j] instanceof SampleDateModel) { final SampleDateModel SAM = ((SampleDateModel) selectedSampleDateModels[i][j]); double wMean = SAM.getValue().movePointLeft(6).doubleValue(); double wMeanOneSigma = SAM.getOneSigmaAbs().movePointLeft(6).doubleValue(); Path2D mean = new Path2D.Double(Path2D.WIND_NON_ZERO); // july 2008 // modified to show de-selected fractions as gray // this means a new special list of fractionIDs is created fromall non-rejected fractions // and each instance is tested for being included // should eventually refactor Vector<String> allFIDs = new Vector<String>(); for (String f : ((UPbReduxAliquot) SAM.getAliquot()).getAliquotFractionIDs()) { // test added for Sample-based wm if (SAM.fractionDateIsPositive(// ((UPbReduxAliquot) SAM.getAliquot()).getAliquotFractionByName(f))) { allFIDs.add(f); } } final int iFinal = i; if (fractionSortOrder.equalsIgnoreCase("weight")) { Collections.sort(allFIDs, new Comparator<String>() { public int compare(String fID1, String fID2) { double invertOneSigmaF1 = // 1.0 // / ((UPbReduxAliquot) selectedSampleDateModels[iFinal][0])// .getAliquotFractionByName(fID1)// .getRadiogenicIsotopeDateByName(SAM.getDateName())// .getOneSigmaAbs().movePointLeft(6).doubleValue(); double invertOneSigmaF2 = // 1.0 // / ((UPbReduxAliquot) selectedSampleDateModels[iFinal][0])// .getAliquotFractionByName(fID2)// .getRadiogenicIsotopeDateByName(SAM.getDateName())// .getOneSigmaAbs().movePointLeft(6).doubleValue(); return Double.compare(invertOneSigmaF2, invertOneSigmaF1); } }); } else if (fractionSortOrder.equalsIgnoreCase("date")) { Collections.sort(allFIDs, new Comparator<String>() { public int compare(String fID1, String fID2) { double dateF1 = // ((UPbReduxAliquot) selectedSampleDateModels[iFinal][0])// .getAliquotFractionByName(fID1)// .getRadiogenicIsotopeDateByName(SAM.getDateName())// .getValue().doubleValue(); double dateF2 = // ((UPbReduxAliquot) selectedSampleDateModels[iFinal][0])// .getAliquotFractionByName(fID2)// .getRadiogenicIsotopeDateByName(SAM.getDateName())// .getValue().doubleValue(); return Double.compare(dateF1, dateF2); } }); } else if ( /* ! isInRandomMode() &&*/fractionSortOrder.equalsIgnoreCase("random")) { Collections.shuffle(allFIDs, new Random()); } else if (fractionSortOrder.equalsIgnoreCase("name")) { // default to alphabetic by name //Collections.sort(allFIDs); // april 2010 give same lexigraphic ordering that UPbFractions get Collections.sort(allFIDs, new IntuitiveStringComparator<String>()); } else { // do nothing } double actualWidthX = (allFIDs.size()) * (barWidth + barGap);//; + barGap; // plot 2-sigma of mean mean.moveTo((float) mapX(startSamX, getMinX_Display(), rangeX, graphWidth), (float) mapY(wMean + 2.0 * wMeanOneSigma, getMaxY_Display(), rangeY, graphHeight)); mean.lineTo((float) mapX(startSamX + actualWidthX, getMinX_Display(), rangeX, graphWidth), (float) mapY(wMean + 2.0 * wMeanOneSigma, getMaxY_Display(), rangeY, graphHeight)); mean.lineTo((float) mapX(startSamX + actualWidthX, getMinX_Display(), rangeX, graphWidth), (float) mapY(wMean - 2.0 * wMeanOneSigma, getMaxY_Display(), rangeY, graphHeight)); mean.lineTo((float) mapX(startSamX, getMinX_Display(), rangeX, graphWidth), (float) mapY(wMean - 2.0 * wMeanOneSigma, getMaxY_Display(), rangeY, graphHeight)); mean.closePath(); g2d.setColor(ReduxConstants.mySampleYellowColor); g2d.fill(mean); g2d.setPaint(Color.BLACK); // plot 1-sigma of mean mean.reset(); mean.moveTo((float) mapX(startSamX, getMinX_Display(), rangeX, graphWidth), (float) mapY(wMean + wMeanOneSigma, getMaxY_Display(), rangeY, graphHeight)); mean.lineTo((float) mapX(startSamX + actualWidthX, getMinX_Display(), rangeX, graphWidth), (float) mapY(wMean + wMeanOneSigma, getMaxY_Display(), rangeY, graphHeight)); mean.lineTo((float) mapX(startSamX + actualWidthX, getMinX_Display(), rangeX, graphWidth), (float) mapY(wMean - wMeanOneSigma, getMaxY_Display(), rangeY, graphHeight)); mean.lineTo((float) mapX(startSamX, getMinX_Display(), rangeX, graphWidth), (float) mapY(wMean - wMeanOneSigma, getMaxY_Display(), rangeY, graphHeight)); mean.closePath(); g2d.setColor(ReduxConstants.ColorOfRedux); g2d.fill(mean); g2d.setPaint(Color.BLACK); // plot mean mean.reset(); mean.moveTo((float) mapX(startSamX, getMinX_Display(), rangeX, graphWidth), (float) mapY(wMean, getMaxY_Display(), rangeY, graphHeight)); mean.lineTo((float) mapX(startSamX + actualWidthX, getMinX_Display(), rangeX, graphWidth), (float) mapY(wMean, getMaxY_Display(), rangeY, graphHeight)); g2d.setStroke(new BasicStroke(1.0f)); g2d.draw(mean); g2d.setStroke(new BasicStroke(2.0f)); saveStartSamX = startSamX; // plot fraction bars double minPoint = 5000.0; double maxWeight = 0.0; double totalWeight = 0.0; int barNum = 0; for (String fID : allFIDs) { // the dateModel has an associated aliquot, but in sample mode, it is a // standin aliquot for the sample. to get the aliquot number for // use in coloring fractions, we need to query the fraction itself String aliquotName = sample.getAliquotNameByFractionID(fID); Color includedFillColor = new Color(0, 0, 0); if (sample.getSampleDateInterpretationGUISettings().getAliquotOptions().get(aliquotName) .containsKey("includedFillColor")) { String[] temp = // sample.getSampleDateInterpretationGUISettings().getAliquotOptions() .get(aliquotName).get("includedFillColor").split(","); includedFillColor = buildRGBColor(temp); } Fraction f = ((UPbReduxAliquot) selectedSampleDateModels[i][0]) .getAliquotFractionByName(fID); double date = f.//((UPbReduxAliquot) selectedSampleDateModels[i][0]).getAliquotFractionByName(fID).// getRadiogenicIsotopeDateByName(SAM.getDateName()).getValue().movePointLeft(6) .doubleValue(); double twoSigma = f.//((UPbReduxAliquot) selectedSampleDateModels[i][0]).getAliquotFractionByName(fID).// getRadiogenicIsotopeDateByName(SAM.getDateName()).getTwoSigmaAbs().movePointLeft(6) .doubleValue(); if ((date - twoSigma) < minPoint) { minPoint = (date - twoSigma); } double invertedOneSigma = // 1.0 // / f.//((UPbReduxAliquot) selectedSampleDateModels[i][0]).getAliquotFractionByName(fID).// getRadiogenicIsotopeDateByName(SAM.getDateName()).getOneSigmaAbs() .movePointLeft(6).doubleValue(); if (invertedOneSigma > maxWeight) { maxWeight = invertedOneSigma; } Path2D bar = new Path2D.Double(Path2D.WIND_NON_ZERO); bar.moveTo( (float) mapX(saveStartSamX + ((barGap / 2.0) + barNum * (barWidth + barGap)), getMinX_Display(), rangeX, graphWidth), (float) mapY(date + twoSigma, getMaxY_Display(), rangeY, graphHeight)); bar.lineTo( (float) mapX( saveStartSamX + ((barGap / 2.0) + barNum * (barWidth + barGap)) + barWidth, getMinX_Display(), rangeX, graphWidth), (float) mapY(date + twoSigma, getMaxY_Display(), rangeY, graphHeight)); bar.lineTo( (float) mapX( saveStartSamX + ((barGap / 2.0) + barNum * (barWidth + barGap)) + barWidth, getMinX_Display(), rangeX, graphWidth), (float) mapY(date - twoSigma, getMaxY_Display(), rangeY, graphHeight)); bar.lineTo( (float) mapX(saveStartSamX + ((barGap / 2.0) + barNum * (barWidth + barGap)), getMinX_Display(), rangeX, graphWidth), (float) mapY(date - twoSigma, getMaxY_Display(), rangeY, graphHeight)); bar.closePath(); Composite originalComposite = g2d.getComposite(); if (SAM.getIncludedFractionIDsVector().contains(fID)) { g2d.setComposite(AlphaComposite.getInstance(AlphaComposite.SRC_OVER, 0.8f)); totalWeight += Math.pow(invertedOneSigma, 2.0); // april 2014 experiment if (f.getRgbColor() != 0) { includedFillColor = new Color(f.getRgbColor()); } } else { g2d.setComposite(AlphaComposite.getInstance(AlphaComposite.SRC_OVER, 0.2f)); } g2d.setPaint(includedFillColor); g2d.draw(bar); //restore composite g2d.setComposite(originalComposite); g2d.setColor(Color.black); // label fraction at top g2d.rotate(-Math.PI / 4.0, (float) mapX(saveStartSamX + ((barGap / 2.0) + barNum * (barWidth + barGap)), getMinX_Display(), rangeX, graphWidth), (float) mapY(date + twoSigma, getMaxY_Display(), rangeY, graphHeight)); g2d.drawString( ((UPbReduxAliquot) selectedSampleDateModels[i][0]).getAliquotFractionByName(fID) .getFractionID(), (float) mapX(saveStartSamX + ((barGap / 2.0) + barNum * (barWidth + barGap)), getMinX_Display(), rangeX, graphWidth) + 15, (float) mapY(date + twoSigma, getMaxY_Display(), rangeY, graphHeight)); g2d.rotate(Math.PI / 4.0, (float) mapX(saveStartSamX + ((barGap / 2.0) + barNum * (barWidth + barGap)), getMinX_Display(), rangeX, graphWidth), (float) mapY(date + twoSigma, getMaxY_Display(), rangeY, graphHeight)); barNum++; // startSamX += 2 * barWidth; startSamX += barWidth + barGap; } // display three info boxes below weighted means // each tic is the height of one calculated y-axis tic // determine the y axis tic double minYtic = Math.ceil(getMinY_Display() * 100) / 100; double maxYtic = Math.floor(getMaxY_Display() * 100) / 100; double deltay = Math.rint((maxYtic - minYtic) * 10 + 0.5); double yTic = deltay / 100; double yTopSummary = minPoint - yTic / 2.0;// wMeanOneSigma; //double specialYTic = yTic; double yTopWeights = yTopSummary - yTic * 1.1; double yTopMSWD_PDF = yTopWeights - yTic * 1.1; // summary box Path2D box = new Path2D.Double(Path2D.WIND_NON_ZERO); box.moveTo((float) mapX(saveStartSamX, getMinX_Display(), rangeX, graphWidth), (float) mapY(yTopSummary, getMaxY_Display(), rangeY, graphHeight)); box.lineTo((float) mapX(saveStartSamX + actualWidthX, getMinX_Display(), rangeX, graphWidth), (float) mapY(yTopSummary, getMaxY_Display(), rangeY, graphHeight)); box.lineTo((float) mapX(saveStartSamX + actualWidthX, getMinX_Display(), rangeX, graphWidth), (float) mapY(yTopSummary - yTic, getMaxY_Display(), rangeY, graphHeight)); box.lineTo((float) mapX(saveStartSamX, getMinX_Display(), rangeX, graphWidth), (float) mapY(yTopSummary - yTic, getMaxY_Display(), rangeY, graphHeight)); box.closePath(); g2d.setStroke(new BasicStroke(1.5f)); g2d.draw(box); // Info Box g2d.drawString(// SAM.getAliquot().getAliquotName(), (float) mapX(saveStartSamX, getMinX_Display(), rangeX, graphWidth) + 4f, (float) mapY(yTopSummary, getMaxY_Display(), rangeY, graphHeight) + 13f); g2d.drawString(// SAM.getName(), (float) mapX(saveStartSamX, getMinX_Display(), rangeX, graphWidth) + 4f, (float) mapY(yTopSummary, getMaxY_Display(), rangeY, graphHeight) + 25f); g2d.drawString(// SAM.FormatValueAndTwoSigmaABSThreeWaysForPublication(6, 2), (float) mapX(saveStartSamX, getMinX_Display(), rangeX, graphWidth) + 4f, (float) mapY(yTopSummary, getMaxY_Display(), rangeY, graphHeight) + 36f); g2d.drawString(// SAM.ShowCustomMSWDwithN(), (float) mapX(saveStartSamX, getMinX_Display(), rangeX, graphWidth) + 4f, (float) mapY(yTopSummary, getMaxY_Display(), rangeY, graphHeight) + 48f); // weights box box.reset(); box.moveTo((float) mapX(saveStartSamX, getMinX_Display(), rangeX, graphWidth), (float) mapY(yTopWeights, getMaxY_Display(), rangeY, graphHeight)); box.lineTo((float) mapX(saveStartSamX + actualWidthX, getMinX_Display(), rangeX, graphWidth), (float) mapY(yTopWeights, getMaxY_Display(), rangeY, graphHeight)); box.lineTo((float) mapX(saveStartSamX + actualWidthX, getMinX_Display(), rangeX, graphWidth), (float) mapY(yTopWeights - yTic, getMaxY_Display(), rangeY, graphHeight)); box.lineTo((float) mapX(saveStartSamX, getMinX_Display(), rangeX, graphWidth), (float) mapY(yTopWeights - yTic, getMaxY_Display(), rangeY, graphHeight)); box.closePath(); g2d.setStroke(new BasicStroke(1.5f)); g2d.draw(box); // plot fraction weights double artificialXRange = allFIDs.size(); double count = 0; //double weightWidth = Math.min(3.0 * barWidth, (yTic / rangeY * graphHeight)) - 15;//yTic;//barWidth * 2.0; double weightWidth = (barWidth + barGap) * 0.9; for (String fID : allFIDs) { // the dateModel has an associated aliquot, but in sample mode, it is a // standin aliquot for the sample. to get the aliquot number for // use in coloring fractions, we need to query the fraction itself String aliquotName = sample.getAliquotNameByFractionID(fID); Fraction f = ((UPbReduxAliquot) selectedSampleDateModels[i][0]) .getAliquotFractionByName(fID); Color includedFillColor = new Color(0, 0, 0); if (sample.getSampleDateInterpretationGUISettings().getAliquotOptions().get(aliquotName) .containsKey("includedFillColor")) { String[] temp = // sample.getSampleDateInterpretationGUISettings().getAliquotOptions() .get(aliquotName).get("includedFillColor").split(","); includedFillColor = buildRGBColor(temp); } double invertOneSigma = // 1.0 // / ((UPbReduxAliquot) selectedSampleDateModels[i][0]) .getAliquotFractionByName(fID)// .getRadiogenicIsotopeDateByName(SAM.getDateName()).getOneSigmaAbs() .movePointLeft(6).doubleValue(); Path2D weight = new Path2D.Double(Path2D.WIND_NON_ZERO); weight.moveTo( (float) mapX(saveStartSamX + (count + 0.5) / artificialXRange * actualWidthX, getMinX_Display(), rangeX, graphWidth) // - (float) (invertOneSigma / maxWeight / 2.0 * weightWidth), (float) mapY(yTopWeights - (yTic / 2.0), getMaxY_Display(), rangeY, graphHeight) // + (float) (invertOneSigma / maxWeight / 2.0 * weightWidth) - 5f); weight.lineTo( (float) mapX(saveStartSamX + (count + 0.5) / artificialXRange * actualWidthX, getMinX_Display(), rangeX, graphWidth) // + (float) (invertOneSigma / maxWeight / 2.0 * weightWidth), (float) mapY(yTopWeights - (yTic / 2.0), getMaxY_Display(), rangeY, graphHeight) // + (float) (invertOneSigma / maxWeight / 2.0 * weightWidth) - 5f); weight.lineTo( (float) mapX(saveStartSamX + (count + 0.5) / artificialXRange * actualWidthX, getMinX_Display(), rangeX, graphWidth) // + (float) (invertOneSigma / maxWeight / 2.0 * weightWidth), (float) mapY(yTopWeights - (yTic / 2.0), getMaxY_Display(), rangeY, graphHeight) // - (float) (invertOneSigma / maxWeight / 2.0 * weightWidth) - 5f); weight.lineTo( (float) mapX(saveStartSamX + (count + 0.5) / artificialXRange * actualWidthX, getMinX_Display(), rangeX, graphWidth) // - (float) (invertOneSigma / maxWeight / 2.0 * weightWidth), (float) mapY(yTopWeights - (yTic / 2.0), getMaxY_Display(), rangeY, graphHeight) // - (float) (invertOneSigma / maxWeight / 2.0 * weightWidth) - 5f); weight.closePath(); g2d.setStroke(new BasicStroke(2.5f)); // test for included or not == black or gray String weightPerCent = " 0";//0.0%"; // g2d.setPaint(includedFillColor); Composite originalComposite = g2d.getComposite(); if (SAM.getIncludedFractionIDsVector().contains(fID)) { g2d.setComposite(AlphaComposite.getInstance(AlphaComposite.SRC_OVER, 0.8f)); weightPerCent = formatter1DecPlace .format(Math.pow(invertOneSigma, 2.0) / totalWeight * 100.0);// + "%"; // april 2014 experiment if (f.getRgbColor() != 0) { includedFillColor = new Color(f.getRgbColor()); } } else { g2d.setComposite(AlphaComposite.getInstance(AlphaComposite.SRC_OVER, 0.2f)); } g2d.setPaint(includedFillColor); g2d.fill(weight); //restore composite g2d.setComposite(originalComposite); // write percent of total weight g2d.drawString(weightPerCent, (float) mapX(saveStartSamX + (count + 0.5) / artificialXRange * actualWidthX, getMinX_Display(), rangeX, graphWidth) // - (float) (invertOneSigma / maxWeight / 2.0 * weightWidth), (float) mapY(yTopWeights - yTic, getMaxY_Display(), rangeY, graphHeight) - 5f); g2d.setColor(Color.black); count += 1.0; } // double box height for graph yTic *= 2.0; // plot MSWD_PDF // store function x,y values Vector<Double> xVals = new Vector<Double>(); Vector<Double> yVals = new Vector<Double>(); double f = SAM.getIncludedFractionIDsVector().size() - 1; if (f > 1.0) { g2d.setStroke(new BasicStroke(1.0f)); double yRange = MSWDCoordinates.valuesByPointCount[(int) f][5] * 1.03; // alitle air at the top of curve double xStart = MSWDCoordinates.valuesByPointCount[(int) f][1]; double xRange = MSWDCoordinates.valuesByPointCount[(int) f][4] - xStart; double xStep = 0.005; Path2D MSWD_PDF = new Path2D.Double(Path2D.WIND_NON_ZERO); Path2D MSWD_right = new Path2D.Double(Path2D.WIND_NON_ZERO); // start at lower left corner of box (may or may not be 0,0 ) MSWD_PDF.moveTo(// (float) mapX((double) saveStartSamX, getMinX_Display(), rangeX, graphWidth), (float) mapY(yTopMSWD_PDF - yTic, getMaxY_Display(), rangeY, graphHeight)); // setup MSWD to paint last Path2D MSWD = null; // calculate function values for (double x = xStart; x < xRange; x += xStep) { xVals.add((((x - xStart) / xRange) * actualWidthX) + (double) saveStartSamX); double y = // Math.pow(2, -1.0 * f / 2.0)// * Math.exp(-1.0 * f * x / 2.0)// * Math.pow(f, f / 2.0)// * Math.pow(x, (-1.0 + f / 2.0))// / Math.exp(Gamma.logGamma(f / 2.0)); yVals.add(((y / yRange) * yTic) + yTopMSWD_PDF - yTic); MSWD_PDF.lineTo(// (float) mapX(xVals.lastElement(), getMinX_Display(), rangeX, graphWidth), (float) mapY(yVals.lastElement(), getMaxY_Display(), rangeY, graphHeight)); // test for location of left RED zone if ((MSWDCoordinates.valuesByPointCount[(int) f][2] >= x) && (MSWDCoordinates.valuesByPointCount[(int) f][2] < (x + xStep))) { double leftX = MSWDCoordinates.valuesByPointCount[(int) f][2]; xVals.add((((leftX - xStart) / xRange) * actualWidthX) + (double) saveStartSamX); double leftY = // Math.pow(2, -1.0 * f / 2.0)// * Math.exp(-1.0 * f * leftX / 2.0)// * Math.pow(f, f / 2.0)// * Math.pow(leftX, (-1.0 + f / 2.0))// / Math.exp(Gamma.logGamma(f / 2.0)); yVals.add(((leftY / yRange) * yTic) + yTopMSWD_PDF - yTic); MSWD_PDF.lineTo(// (float) mapX(xVals.lastElement(), getMinX_Display(), rangeX, graphWidth), (float) mapY(yVals.lastElement(), getMaxY_Display(), rangeY, graphHeight)); Path2D ciLower = new Path2D.Double(Path2D.WIND_NON_ZERO); ciLower.append(MSWD_PDF.getPathIterator(new AffineTransform()), true); ciLower.lineTo(// (float) mapX(xVals.lastElement(), getMinX_Display(), rangeX, graphWidth), (float) mapY(yTopMSWD_PDF - yTic, getMaxY_Display(), rangeY, graphHeight)); ciLower.closePath(); g2d.setColor(Color.RED); g2d.fill(ciLower); // draw right hand border line to compensate for a bug in the filler Line2D right = new Line2D.Double(// mapX(xVals.lastElement(), getMinX_Display(), rangeX, graphWidth), mapY(yVals.lastElement(), getMaxY_Display(), rangeY, graphHeight), mapX(xVals.lastElement(), getMinX_Display(), rangeX, graphWidth), mapY(yTopMSWD_PDF - yTic, getMaxY_Display(), rangeY, graphHeight)); g2d.setStroke(new BasicStroke(0.5f)); g2d.draw(right); g2d.setStroke(new BasicStroke(1.0f)); g2d.setColor(Color.BLACK); System.out.println("Left Red = (" + leftX + ", " + leftY + ")"); } // test for location of right RED zone if ((MSWDCoordinates.valuesByPointCount[(int) f][3] >= x) && (MSWDCoordinates.valuesByPointCount[(int) f][3] < (x + xStep))) { double rightX = MSWDCoordinates.valuesByPointCount[(int) f][3]; xVals.add((((rightX - xStart) / xRange) * actualWidthX) + (double) saveStartSamX); double rightY = // Math.pow(2, -1.0 * f / 2.0)// * Math.exp(-1.0 * f * rightX / 2.0)// * Math.pow(f, f / 2.0)// * Math.pow(rightX, (-1.0 + f / 2.0))// / Math.exp(Gamma.logGamma(f / 2.0)); yVals.add(((rightY / yRange) * yTic) + yTopMSWD_PDF - yTic); MSWD_PDF.lineTo(// (float) mapX(xVals.lastElement(), getMinX_Display(), rangeX, graphWidth), (float) mapY(yVals.lastElement(), getMaxY_Display(), rangeY, graphHeight)); // here the strategy is to draw the curve and then reset it to record the remainder g2d.setStroke(new BasicStroke(1.0f)); g2d.draw(MSWD_PDF); MSWD_PDF = new Path2D.Double(Path2D.WIND_NON_ZERO); MSWD_PDF.moveTo(// (float) mapX(xVals.lastElement(), getMinX_Display(), rangeX, graphWidth), (float) mapY(yVals.lastElement(), getMaxY_Display(), rangeY, graphHeight)); MSWD_right.moveTo(// (float) mapX(xVals.lastElement(), getMinX_Display(), rangeX, graphWidth), (float) mapY(yTopMSWD_PDF - yTic, getMaxY_Display(), rangeY, graphHeight)); MSWD_right.lineTo(// (float) mapX(xVals.lastElement(), getMinX_Display(), rangeX, graphWidth), (float) mapY(yVals.lastElement(), getMaxY_Display(), rangeY, graphHeight)); System.out.println("Right Red = (" + rightX + ", " + rightY + ")"); } // test for location of MSWD AND paint last if ((SAM.getMeanSquaredWeightedDeviation().doubleValue() >= x) && (SAM.getMeanSquaredWeightedDeviation().doubleValue() < (x + xStep))) { MSWD = new Path2D.Double(Path2D.WIND_NON_ZERO); MSWD.moveTo(// (float) mapX(xVals.lastElement(), getMinX_Display(), rangeX, graphWidth), (float) mapY(yTopMSWD_PDF - yTic, getMaxY_Display(), rangeY, graphHeight)); MSWD.lineTo(// (float) mapX(xVals.lastElement(), getMinX_Display(), rangeX, graphWidth), (float) mapY(yVals.lastElement(), getMaxY_Display(), rangeY, graphHeight)); } } g2d.setStroke(new BasicStroke(1.0f)); // merge with border of right RED and fill MSWD_right.append(MSWD_PDF.getPathIterator(new AffineTransform()), true); g2d.setColor(Color.RED); g2d.fill(MSWD_right); g2d.setColor(Color.BLACK); // draw the remaining curves g2d.draw(MSWD_PDF); // MSWD may be off the graph and hence not exist try { g2d.draw(MSWD); } catch (Exception e) { } // label 95% conf interval and MSWD g2d.drawString(// "95% CI: (" + formatter2DecPlaces.format(MSWDCoordinates.valuesByPointCount[(int) f][2]) + ", " + formatter2DecPlaces.format(MSWDCoordinates.valuesByPointCount[(int) f][3]) + ")", (float) mapX(saveStartSamX + (actualWidthX / 2.0), getMinX_Display(), rangeX, graphWidth) - 30f, (float) mapY(yTopMSWD_PDF, getMaxY_Display(), rangeY, graphHeight) + 15f); // determine if MSWD is out of range String mswdAlert = ""; if (SAM.getMeanSquaredWeightedDeviation() .doubleValue() > MSWDCoordinates.valuesByPointCount[(int) f][4]) { mswdAlert = "\n !Out of Range!"; } g2d.drawString(// "MSWD = " + formatter2DecPlaces .format(SAM.getMeanSquaredWeightedDeviation().doubleValue()) + ", n = " + (int) (f + 1) + mswdAlert, (float) mapX(saveStartSamX + (actualWidthX / 2.0), getMinX_Display(), rangeX, graphWidth) - 15f, (float) mapY(yTopMSWD_PDF, getMaxY_Display(), rangeY, graphHeight) + 30f); } else { g2d.drawString("need more data...", (float) mapX((double) saveStartSamX, getMinX_Display(), rangeX, graphWidth) + 4f, (float) mapY(yTopMSWD_PDF - yTic, getMaxY_Display(), rangeY, graphHeight) - 10f); } // MSWD_PDF box box.reset(); box.moveTo((float) mapX(saveStartSamX, getMinX_Display(), rangeX, graphWidth), (float) mapY(yTopMSWD_PDF, getMaxY_Display(), rangeY, graphHeight)); box.lineTo((float) mapX(saveStartSamX + actualWidthX, getMinX_Display(), rangeX, graphWidth), (float) mapY(yTopMSWD_PDF, getMaxY_Display(), rangeY, graphHeight)); box.lineTo((float) mapX(saveStartSamX + actualWidthX, getMinX_Display(), rangeX, graphWidth), (float) mapY(yTopMSWD_PDF - yTic, getMaxY_Display(), rangeY, graphHeight)); box.lineTo((float) mapX(saveStartSamX, getMinX_Display(), rangeX, graphWidth), (float) mapY(yTopMSWD_PDF - yTic, getMaxY_Display(), rangeY, graphHeight)); box.closePath(); g2d.setStroke(new BasicStroke(1.5f)); g2d.draw(box); // MSWD_PDF x-axis tics if (f > 1.0) { g2d.setStroke(new BasicStroke(1.0f)); double xStart = (MSWDCoordinates.valuesByPointCount[(int) f][1] <= 0.5) ? 0.5 : 1.0; double xRange = MSWDCoordinates.valuesByPointCount[(int) f][4] - MSWDCoordinates.valuesByPointCount[(int) f][1]; double xStep = 0.5; for (double x = xStart; x < xRange; x += xStep) { double xPlot = (((x - MSWDCoordinates.valuesByPointCount[(int) f][1]) / xRange) * actualWidthX) + (double) saveStartSamX; Line2D line = new Line2D.Double(mapX(xPlot, getMinX_Display(), rangeX, graphWidth), mapY(yTopMSWD_PDF - yTic, getMaxY_Display(), rangeY, graphHeight), mapX(xPlot, getMinX_Display(), rangeX, graphWidth), mapY(yTopMSWD_PDF - yTic, getMaxY_Display(), rangeY, graphHeight) + 7); g2d.draw(line); g2d.rotate(-Math.PI / 2.0, (float) mapX(xPlot, getMinX_Display(), rangeX, graphWidth), (float) mapY(yTopMSWD_PDF - yTic, getMaxY_Display(), rangeY, graphHeight)); g2d.drawString(formatter1DecPlace.format(x), (float) mapX(xPlot, getMinX_Display(), rangeX, graphWidth) - 30f, (float) mapY(yTopMSWD_PDF - yTic, getMaxY_Display(), rangeY, graphHeight) + 5f); g2d.rotate(Math.PI / 2.0, (float) mapX(xPlot, getMinX_Display(), rangeX, graphWidth), (float) mapY(yTopMSWD_PDF - yTic, getMaxY_Display(), rangeY, graphHeight)); } } // set counters barNum += samSpace; startSamX += 2 * samSpace * barWidth; } } } // // prevents re-randomization // setInRandomMode( true ); drawAxesAndTicks(g2d, rangeX, rangeY); // draw zoom box if in use if ((Math.abs(zoomMaxX - zoomMinX) * Math.abs(zoomMinY - zoomMaxY)) > 0.0) { g2d.setStroke(new BasicStroke(2.0f)); g2d.setColor(Color.red); g2d.drawRect(// Math.min(zoomMinX, zoomMaxX), Math.min(zoomMaxY, zoomMinY), Math.abs(zoomMaxX - zoomMinX), Math.abs(zoomMinY - zoomMaxY)); } }
From source file:org.genemania.engine.actions.ComputeEnrichment.java
public static double computeHyperGeo(double x, double N, double n, double k) { double h = Gamma.logGamma(k + 1) - Gamma.logGamma(k - x + 1) - Gamma.logGamma(x + 1) + Gamma.logGamma(N - k + 1) - Gamma.logGamma(N - k - n + x + 1) - Gamma.logGamma(n - x + 1) - Gamma.logGamma(N + 1) + Gamma.logGamma(N - n + 1) + Gamma.logGamma(n + 1); return Math.exp(h); }