Examples with GammaDistribution org.apache.commons.math3.distribution.GammaDistribution used on opensource projects

Example 1 with GammaDistribution

use of org.apache.commons.math3.distribution.GammaDistribution in project GDSC-SMLM by aherbert.

the class CreateData method createPhotonDistribution.

/**
	 * @return A photon distribution loaded from a file of floating-point values with the specified population mean.
	 */
private RealDistribution createPhotonDistribution() {
    if (PHOTON_DISTRIBUTION[PHOTON_CUSTOM].equals(settings.photonDistribution)) {
        // Get the distribution file
        String filename = Utils.getFilename("Photon_distribution", settings.photonDistributionFile);
        if (filename != null) {
            settings.photonDistributionFile = filename;
            try {
                InputStream is = new FileInputStream(new File(settings.photonDistributionFile));
                BufferedReader in = new BufferedReader(new UnicodeReader(is, null));
                StoredDataStatistics stats = new StoredDataStatistics();
                try {
                    String str = null;
                    double val = 0.0d;
                    while ((str = in.readLine()) != null) {
                        val = Double.parseDouble(str);
                        stats.add(val);
                    }
                } finally {
                    in.close();
                }
                if (stats.getSum() > 0) {
                    // Update the statistics to the desired mean.
                    double scale = (double) settings.photonsPerSecond / stats.getMean();
                    double[] values = stats.getValues();
                    for (int i = 0; i < values.length; i++) values[i] *= scale;
                    // TODO - Investigate the limits of this distribution. 
                    // How far above and below the input data will values be generated.
                    // Create the distribution using the recommended number of bins
                    final int binCount = stats.getN() / 10;
                    EmpiricalDistribution dist = new EmpiricalDistribution(binCount, createRandomGenerator());
                    dist.load(values);
                    return dist;
                }
            } catch (IOException e) {
            // Ignore
            } catch (NullArgumentException e) {
            // Ignore 
            } catch (NumberFormatException e) {
            // Ignore
            }
        }
        Utils.log("Failed to load custom photon distribution from file: %s. Default to fixed.", settings.photonDistributionFile);
    } else if (PHOTON_DISTRIBUTION[PHOTON_UNIFORM].equals(settings.photonDistribution)) {
        if (settings.photonsPerSecond < settings.photonsPerSecondMaximum) {
            UniformRealDistribution dist = new UniformRealDistribution(createRandomGenerator(), settings.photonsPerSecond, settings.photonsPerSecondMaximum);
            return dist;
        }
    } else if (PHOTON_DISTRIBUTION[PHOTON_GAMMA].equals(settings.photonDistribution)) {
        final double scaleParameter = settings.photonsPerSecond / settings.photonShape;
        GammaDistribution dist = new GammaDistribution(createRandomGenerator(), settings.photonShape, scaleParameter, ExponentialDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
        return dist;
    } else if (PHOTON_DISTRIBUTION[PHOTON_CORRELATED].equals(settings.photonDistribution)) {
        // No distribution required
        return null;
    }
    settings.photonDistribution = PHOTON_DISTRIBUTION[PHOTON_FIXED];
    return null;
}

Example 2 with GammaDistribution

use of org.apache.commons.math3.distribution.GammaDistribution in project gatk by broadinstitute.

the class AlleleFractionSegmenterUnitTest method generateCounts.

//visible for testing joint segmentation
protected static AllelicCountCollection generateCounts(final List<Double> minorAlleleFractionSequence, final List<SimpleInterval> positions, final RandomGenerator rng, final AlleleFractionGlobalParameters trueParams) {
    //translate to ApacheCommons' parametrization of the gamma distribution
    final GammaDistribution biasGenerator = getGammaDistribution(trueParams, rng);
    final double outlierProbability = trueParams.getOutlierProbability();
    final AllelicCountCollection counts = new AllelicCountCollection();
    for (int n = 0; n < minorAlleleFractionSequence.size(); n++) {
        counts.add(generateAllelicCount(minorAlleleFractionSequence.get(n), positions.get(n), rng, biasGenerator, outlierProbability));
    }
    return counts;
}

Example 3 with GammaDistribution

use of org.apache.commons.math3.distribution.GammaDistribution in project gatk by broadinstitute.

the class AlleleFractionSegmenterUnitTest method generateAllelicCount.

protected static AllelicCount generateAllelicCount(final double minorFraction, final SimpleInterval position, final RandomGenerator rng, final GammaDistribution biasGenerator, final double outlierProbability) {
    final int numReads = 100;
    final double bias = biasGenerator.sample();
    //flip a coin to decide alt minor (alt fraction = minor fraction) or ref minor (alt fraction = 1 - minor fraction)
    final double altFraction = rng.nextDouble() < 0.5 ? minorFraction : 1 - minorFraction;
    //the probability of an alt read is the alt fraction modified by the bias or, in the case of an outlier, random
    final double pAlt = rng.nextDouble() < outlierProbability ? rng.nextDouble() : altFraction / (altFraction + (1 - altFraction) * bias);
    final int numAltReads = new BinomialDistribution(rng, numReads, pAlt).sample();
    final int numRefReads = numReads - numAltReads;
    return new AllelicCount(position, numAltReads, numRefReads);
}

Example 4 with GammaDistribution

use of org.apache.commons.math3.distribution.GammaDistribution in project gatk-protected by broadinstitute.

the class AlleleFractionSegmenterUnitTest method generateCounts.

//visible for testing joint segmentation
protected static AllelicCountCollection generateCounts(final List<Double> minorAlleleFractionSequence, final List<SimpleInterval> positions, final RandomGenerator rng, final AlleleFractionGlobalParameters trueParams) {
    //translate to ApacheCommons' parametrization of the gamma distribution
    final GammaDistribution biasGenerator = getGammaDistribution(trueParams, rng);
    final double outlierProbability = trueParams.getOutlierProbability();
    final AllelicCountCollection counts = new AllelicCountCollection();
    for (int n = 0; n < minorAlleleFractionSequence.size(); n++) {
        counts.add(generateAllelicCount(minorAlleleFractionSequence.get(n), positions.get(n), rng, biasGenerator, outlierProbability));
    }
    return counts;
}

Example 5 with GammaDistribution

use of org.apache.commons.math3.distribution.GammaDistribution in project GDSC-SMLM by aherbert.

the class DiffusionRateTest method plotJumpDistances.

/**
	 * Plot a cumulative histogram and standard histogram of the jump distances.
	 *
	 * @param title
	 *            the title
	 * @param jumpDistances
	 *            the jump distances
	 * @param dimensions
	 *            the number of dimensions for the jumps
	 * @param steps
	 *            the steps
	 */
private void plotJumpDistances(String title, DoubleData jumpDistances, int dimensions) {
    // Cumulative histogram
    // --------------------
    double[] values = jumpDistances.values();
    String title2 = title + " Cumulative Jump Distance " + dimensions + "D";
    double[][] jdHistogram = JumpDistanceAnalysis.cumulativeHistogram(values);
    Plot2 jdPlot = new Plot2(title2, "Distance (um^2)", "Cumulative Probability", jdHistogram[0], jdHistogram[1]);
    PlotWindow pw2 = Utils.display(title2, jdPlot);
    if (Utils.isNewWindow())
        idList[idCount++] = pw2.getImagePlus().getID();
    // Plot the expected function
    // This is the Chi-squared distribution: The sum of the squares of k independent
    // standard normal random variables with k = dimensions. It is a special case of
    // the gamma distribution. If the normals have non-unit variance the distribution 
    // is scaled.
    // Chi       ~ Gamma(k/2, 2)      // using the scale parameterisation of the gamma
    // s^2 * Chi ~ Gamma(k/2, 2*s^2)
    // So if s^2 = 2D:
    // 2D * Chi  ~ Gamma(k/2, 4D)
    double estimatedD = simpleD * simpleSteps;
    double max = Maths.max(values);
    double[] x = Utils.newArray(1000, 0, max / 1000);
    double k = dimensions / 2.0;
    double mean = 4 * estimatedD;
    GammaDistribution dist = new GammaDistribution(k, mean);
    double[] y = new double[x.length];
    for (int i = 0; i < x.length; i++) y[i] = dist.cumulativeProbability(x[i]);
    jdPlot.setColor(Color.red);
    jdPlot.addPoints(x, y, Plot.LINE);
    Utils.display(title2, jdPlot);
    // Histogram
    // ---------
    title2 = title + " Jump " + dimensions + "D";
    int plotId = Utils.showHistogram(title2, jumpDistances, "Distance (um^2)", 0, 0, Math.max(20, values.length / 1000));
    if (Utils.isNewWindow())
        idList[idCount++] = plotId;
    // Recompute the expected function
    for (int i = 0; i < x.length; i++) y[i] = dist.density(x[i]);
    // Scale to have the same area
    if (Utils.xValues.length > 1) {
        final double area1 = jumpDistances.size() * (Utils.xValues[1] - Utils.xValues[0]);
        final double area2 = dist.cumulativeProbability(x[x.length - 1]);
        final double scaleFactor = area1 / area2;
        for (int i = 0; i < y.length; i++) y[i] *= scaleFactor;
    }
    jdPlot = Utils.plot;
    jdPlot.setColor(Color.red);
    jdPlot.addPoints(x, y, Plot.LINE);
    Utils.display(WindowManager.getImage(plotId).getTitle(), jdPlot);
}