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

Example 1 with BinomialDistribution

use of org.apache.commons.math3.distribution.BinomialDistribution in project deeplearning4j by deeplearning4j.

the class TestReconstructionDistributions method testBernoulliLogProb.

@Test
public void testBernoulliLogProb() {
    Nd4j.getRandom().setSeed(12345);
    int inputSize = 4;
    int[] mbs = new int[] { 1, 2, 5 };
    Random r = new Random(12345);
    for (boolean average : new boolean[] { true, false }) {
        for (int minibatch : mbs) {
            INDArray x = Nd4j.zeros(minibatch, inputSize);
            for (int i = 0; i < minibatch; i++) {
                for (int j = 0; j < inputSize; j++) {
                    x.putScalar(i, j, r.nextInt(2));
                }
            }
            //i.e., pre-sigmoid prob
            INDArray distributionParams = Nd4j.rand(minibatch, inputSize).muli(2).subi(1);
            INDArray prob = Transforms.sigmoid(distributionParams, true);
            ReconstructionDistribution dist = new BernoulliReconstructionDistribution("sigmoid");
            double negLogProb = dist.negLogProbability(x, distributionParams, average);
            INDArray exampleNegLogProb = dist.exampleNegLogProbability(x, distributionParams);
            assertArrayEquals(new int[] { minibatch, 1 }, exampleNegLogProb.shape());
            //Calculate the same thing, but using Apache Commons math
            double logProbSum = 0.0;
            for (int i = 0; i < minibatch; i++) {
                double exampleSum = 0.0;
                for (int j = 0; j < inputSize; j++) {
                    double p = prob.getDouble(i, j);
                    //Bernoulli is a special case of binomial
                    BinomialDistribution binomial = new BinomialDistribution(1, p);
                    double xVal = x.getDouble(i, j);
                    double thisLogProb = binomial.logProbability((int) xVal);
                    logProbSum += thisLogProb;
                    exampleSum += thisLogProb;
                }
                assertEquals(-exampleNegLogProb.getDouble(i), exampleSum, 1e-6);
            }
            double expNegLogProb;
            if (average) {
                expNegLogProb = -logProbSum / minibatch;
            } else {
                expNegLogProb = -logProbSum;
            }
            //                System.out.println(x);
            //                System.out.println(expNegLogProb + "\t" + logProb + "\t" + (logProb / expNegLogProb));
            assertEquals(expNegLogProb, negLogProb, 1e-6);
            //Also: check random sampling...
            int count = minibatch * inputSize;
            INDArray arr = Nd4j.linspace(-3, 3, count).reshape(minibatch, inputSize);
            INDArray sampleMean = dist.generateAtMean(arr);
            INDArray sampleRandom = dist.generateRandom(arr);
            for (int i = 0; i < minibatch; i++) {
                for (int j = 0; j < inputSize; j++) {
                    double d1 = sampleMean.getDouble(i, j);
                    double d2 = sampleRandom.getDouble(i, j);
                    //Mean value - probability... could do 0 or 1 (based on most likely) but that isn't very useful...
                    assertTrue(d1 >= 0.0 || d1 <= 1.0);
                    assertTrue(d2 == 0.0 || d2 == 1.0);
                }
            }
        }
    }
}

Example 2 with BinomialDistribution

use of org.apache.commons.math3.distribution.BinomialDistribution in project presto by prestodb.

the class TestStochasticPriorityQueue method testPollDistribution.

@Test
public void testPollDistribution() {
    StochasticPriorityQueue<String> queue = new StochasticPriorityQueue<>();
    for (int i = 0; i < 100; i++) {
        assertTrue(queue.addOrUpdate("foo" + i, 1));
    }
    for (int i = 0; i < 100; i++) {
        assertTrue(queue.addOrUpdate("bar" + i, 1));
    }
    int foo = 0;
    for (int i = 0; i < 1000; i++) {
        String value = queue.poll();
        if (value.startsWith("foo")) {
            foo++;
        }
        assertTrue(queue.addOrUpdate(value, 1));
    }
    BinomialDistribution binomial = new BinomialDistribution(1000, 0.5);
    int lowerBound = binomial.inverseCumulativeProbability(0.000001);
    int upperBound = binomial.inverseCumulativeProbability(0.999999);
    assertLessThan(foo, upperBound);
    assertGreaterThan(foo, lowerBound);
    // Update foo weights to 2:1 distribution
    for (int i = 0; i < 100; i++) {
        assertFalse(queue.addOrUpdate("foo" + i, 2));
    }
    foo = 0;
    for (int i = 0; i < 1000; i++) {
        String value = queue.poll();
        if (value.startsWith("foo")) {
            foo++;
            assertTrue(queue.addOrUpdate(value, 2));
        } else {
            assertTrue(queue.addOrUpdate(value, 1));
        }
    }
    binomial = new BinomialDistribution(1000, 2.0 / 3.0);
    lowerBound = binomial.inverseCumulativeProbability(0.000001);
    upperBound = binomial.inverseCumulativeProbability(0.999999);
    assertLessThan(foo, upperBound);
    assertGreaterThan(foo, lowerBound);
}

Example 3 with BinomialDistribution

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

the class PulseActivationAnalysis method simulateActivations.

private void simulateActivations(RandomDataGenerator rdg, float[][][] molecules, int c, MemoryPeakResults[] results) {
    int n = molecules[c].length;
    if (n == 0)
        return;
    // Compute desired number per frame
    double umPerPixel = sim_nmPerPixel / 1000;
    double um2PerPixel = umPerPixel * umPerPixel;
    double area = sim_size * sim_size * um2PerPixel;
    double nPerFrame = area * sim_activationDensity;
    // Compute the activation probability (but set an upper limit so not all are on in every frame)
    double p = Math.min(0.5, nPerFrame / n);
    // Determine the other channels activation probability using crosstalk
    double[] p0 = { p, p, p };
    int index1, index2, c1, c2;
    switch(c) {
        case 0:
            index1 = C12;
            index2 = C13;
            c1 = 1;
            c2 = 2;
            break;
        case 1:
            index1 = C21;
            index2 = C23;
            c1 = 0;
            c2 = 2;
            break;
        case 2:
        default:
            index1 = C31;
            index2 = C32;
            c1 = 0;
            c2 = 1;
            break;
    }
    p0[c1] *= ct[index1];
    p0[c2] *= ct[index2];
    // Assume 10 frames after each channel pulse => 30 frames per cycle
    double precision = sim_precision[c] / sim_nmPerPixel;
    int id = c + 1;
    RandomGenerator rand = rdg.getRandomGenerator();
    BinomialDistribution[] bd = new BinomialDistribution[4];
    for (int i = 0; i < 3; i++) bd[i] = createBinomialDistribution(rand, n, p0[i]);
    int[] frames = new int[27];
    for (int i = 1, j = 0; i <= 30; i++) {
        if (i % 10 == 1)
            // Skip specific activation frames 
            continue;
        frames[j++] = i;
    }
    bd[3] = createBinomialDistribution(rand, n, p * sim_nonSpecificFrequency);
    // Count the actual cross talk
    int[] count = new int[3];
    for (int i = 0, t = 1; i < sim_cycles; i++, t += 30) {
        count[0] += simulateActivations(rdg, bd[0], molecules[c], results[c], t, precision, id);
        count[1] += simulateActivations(rdg, bd[1], molecules[c], results[c], t + 10, precision, id);
        count[2] += simulateActivations(rdg, bd[2], molecules[c], results[c], t + 20, precision, id);
        // Add non-specific activations
        if (bd[3] != null) {
            for (int t2 : frames) simulateActivations(rdg, bd[3], molecules[c], results[c], t2, precision, id);
        }
    }
    // Report simulated cross talk
    double[] crosstalk = computeCrosstalk(count, c);
    Utils.log("Simulated crosstalk C%s  %s=>%s, C%s  %s=>%s", ctNames[index1], Utils.rounded(ct[index1]), Utils.rounded(crosstalk[c1]), ctNames[index2], Utils.rounded(ct[index2]), Utils.rounded(crosstalk[c2]));
}

Example 4 with BinomialDistribution

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

the class PCPALMClusters method run.

/*
	 * (non-Javadoc)
	 * 
	 * @see ij.plugin.PlugIn#run(java.lang.String)
	 */
public void run(String arg) {
    SMLMUsageTracker.recordPlugin(this.getClass(), arg);
    if (!showDialog())
        return;
    PCPALMMolecules.logSpacer();
    Utils.log(TITLE);
    PCPALMMolecules.logSpacer();
    long start = System.currentTimeMillis();
    HistogramData histogramData;
    if (fileInput) {
        histogramData = loadHistogram(histogramFile);
    } else {
        histogramData = doClustering();
    }
    if (histogramData == null)
        return;
    float[][] hist = histogramData.histogram;
    // Create a histogram of the cluster sizes
    String title = TITLE + " Molecules/cluster";
    String xTitle = "Molecules/cluster";
    String yTitle = "Frequency";
    // Create the data required for fitting and plotting
    float[] xValues = Utils.createHistogramAxis(hist[0]);
    float[] yValues = Utils.createHistogramValues(hist[1]);
    // Plot the histogram
    float yMax = Maths.max(yValues);
    Plot2 plot = new Plot2(title, xTitle, yTitle, xValues, yValues);
    if (xValues.length > 0) {
        double xPadding = 0.05 * (xValues[xValues.length - 1] - xValues[0]);
        plot.setLimits(xValues[0] - xPadding, xValues[xValues.length - 1] + xPadding, 0, yMax * 1.05);
    }
    Utils.display(title, plot);
    HistogramData noiseData = loadNoiseHistogram(histogramData);
    if (noiseData != null) {
        if (subtractNoise(histogramData, noiseData)) {
            // Update the histogram
            title += " (noise subtracted)";
            xValues = Utils.createHistogramAxis(hist[0]);
            yValues = Utils.createHistogramValues(hist[1]);
            yMax = Maths.max(yValues);
            plot = new Plot2(title, xTitle, yTitle, xValues, yValues);
            if (xValues.length > 0) {
                double xPadding = 0.05 * (xValues[xValues.length - 1] - xValues[0]);
                plot.setLimits(xValues[0] - xPadding, xValues[xValues.length - 1] + xPadding, 0, yMax * 1.05);
            }
            Utils.display(title, plot);
            // Automatically save
            if (autoSave) {
                String newFilename = Utils.replaceExtension(histogramData.filename, ".noise.tsv");
                if (saveHistogram(histogramData, newFilename)) {
                    Utils.log("Saved noise-subtracted histogram to " + newFilename);
                }
            }
        }
    }
    // Fit the histogram
    double[] fitParameters = fitBinomial(histogramData);
    if (fitParameters != null) {
        // Add the binomial to the histogram
        int n = (int) fitParameters[0];
        double p = fitParameters[1];
        Utils.log("Optimal fit : N=%d, p=%s", n, Utils.rounded(p));
        BinomialDistribution dist = new BinomialDistribution(n, p);
        // A zero-truncated binomial was fitted.
        // pi is the adjustment factor for the probability density.
        double pi = 1 / (1 - dist.probability(0));
        if (!fileInput) {
            // Calculate the estimated number of clusters from the observed molecules:
            // Actual = (Observed / p-value) / N
            final double actual = (nMolecules / p) / n;
            Utils.log("Estimated number of clusters : (%d / %s) / %d = %s", nMolecules, Utils.rounded(p), n, Utils.rounded(actual));
        }
        double[] x = new double[n + 2];
        double[] y = new double[n + 2];
        // Scale the values to match those on the histogram
        final double normalisingFactor = count * pi;
        for (int i = 0; i <= n; i++) {
            x[i] = i + 0.5;
            y[i] = dist.probability(i) * normalisingFactor;
        }
        x[n + 1] = n + 1.5;
        y[n + 1] = 0;
        // Redraw the plot since the limits may have changed
        plot = new Plot2(title, xTitle, yTitle, xValues, yValues);
        double xPadding = 0.05 * (xValues[xValues.length - 1] - xValues[0]);
        plot.setLimits(xValues[0] - xPadding, xValues[xValues.length - 1] + xPadding, 0, Maths.maxDefault(yMax, y) * 1.05);
        plot.setColor(Color.magenta);
        plot.addPoints(x, y, Plot2.LINE);
        plot.addPoints(x, y, Plot2.CIRCLE);
        plot.setColor(Color.black);
        Utils.display(title, plot);
    }
    double seconds = (System.currentTimeMillis() - start) / 1000.0;
    String msg = TITLE + " complete : " + seconds + "s";
    IJ.showStatus(msg);
    Utils.log(msg);
    return;
}

Example 5 with BinomialDistribution

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

the class ArtifactStatisticsScorer method calculateArtifactPValue.

/** p-value for being an artifact
     *
     * @param totalAltAlleleCount total number of alt reads
     * @param artifactAltAlleleCount alt read counts in a pair orientation that supports an artifact
     * @param biasP believed bias p value for a binomial distribution in artifact variants
     * @return p-value for this being an artifact
     */
public static double calculateArtifactPValue(final int totalAltAlleleCount, final int artifactAltAlleleCount, final double biasP) {
    ParamUtils.isPositiveOrZero(biasP, "bias parameter must be positive or zero.");
    ParamUtils.isPositiveOrZero(totalAltAlleleCount, "total alt allele count must be positive or zero.");
    ParamUtils.isPositiveOrZero(artifactAltAlleleCount, "artifact supporting alt allele count must be positive or zero.");
    ParamUtils.isPositiveOrZero(totalAltAlleleCount - artifactAltAlleleCount, "Total alt count must be same or greater than the artifact alt count.");
    return new BinomialDistribution(null, totalAltAlleleCount, biasP).cumulativeProbability(artifactAltAlleleCount);
}