Search in sources :

Example 1 with Variance

use of org.apache.commons.math3.stat.descriptive.moment.Variance in project presto by prestodb.

the class TestDoubleVariancePopAggregation method getExpectedValue.

@Override
public Number getExpectedValue(int start, int length) {
    if (length == 0) {
        return null;
    }
    double[] values = new double[length];
    for (int i = 0; i < length; i++) {
        values[i] = start + i;
    }
    Variance variance = new Variance(false);
    return variance.evaluate(values);
}
Also used : Variance(org.apache.commons.math3.stat.descriptive.moment.Variance)

Example 2 with Variance

use of org.apache.commons.math3.stat.descriptive.moment.Variance in project presto by prestodb.

the class TestLongVarianceAggregation method getExpectedValue.

@Override
public Number getExpectedValue(int start, int length) {
    if (length < 2) {
        return null;
    }
    double[] values = new double[length];
    for (int i = 0; i < length; i++) {
        values[i] = start + i;
    }
    Variance variance = new Variance();
    return variance.evaluate(values);
}
Also used : Variance(org.apache.commons.math3.stat.descriptive.moment.Variance)

Example 3 with Variance

use of org.apache.commons.math3.stat.descriptive.moment.Variance in project GDSC-SMLM by aherbert.

the class CMOSAnalysis method simulate.

private void simulate() {
    // Create the offset, variance and gain for each pixel
    int n = size * size;
    float[] pixelOffset = new float[n];
    float[] pixelVariance = new float[n];
    float[] pixelGain = new float[n];
    IJ.showStatus("Creating random per-pixel readout");
    long start = System.currentTimeMillis();
    RandomGenerator rg = new Well19937c();
    PoissonDistribution pd = new PoissonDistribution(rg, offset, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);
    ExponentialDistribution ed = new ExponentialDistribution(rg, variance, ExponentialDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
    totalProgress = n;
    stepProgress = Utils.getProgressInterval(totalProgress);
    for (int i = 0; i < n; i++) {
        if (i % n == 0)
            IJ.showProgress(i, n);
        // Q. Should these be clipped to a sensible range?
        pixelOffset[i] = (float) pd.sample();
        pixelVariance[i] = (float) ed.sample();
        pixelGain[i] = (float) (gain + rg.nextGaussian() * gainSD);
    }
    IJ.showProgress(1);
    // Avoid all the file saves from updating the progress bar and status line
    Utils.setShowStatus(false);
    Utils.setShowProgress(false);
    JLabel statusLine = Utils.getStatusLine();
    progressBar = Utils.getProgressBar();
    // Save to the directory as a stack
    ImageStack simulationStack = new ImageStack(size, size);
    simulationStack.addSlice("Offset", pixelOffset);
    simulationStack.addSlice("Variance", pixelVariance);
    simulationStack.addSlice("Gain", pixelGain);
    simulationImp = new ImagePlus("PerPixel", simulationStack);
    // Only the info property is saved to the TIFF file
    simulationImp.setProperty("Info", String.format("Offset=%s; Variance=%s; Gain=%s +/- %s", Utils.rounded(offset), Utils.rounded(variance), Utils.rounded(gain), Utils.rounded(gainSD)));
    IJ.save(simulationImp, new File(directory, "perPixelSimulation.tif").getPath());
    // Create thread pool and workers
    ExecutorService executor = Executors.newFixedThreadPool(getThreads());
    TurboList<Future<?>> futures = new TurboList<Future<?>>(nThreads);
    // Simulate the zero exposure input.
    // Simulate 20 - 200 photon images.
    int[] photons = new int[] { 0, 20, 50, 100, 200 };
    totalProgress = photons.length * frames;
    stepProgress = Utils.getProgressInterval(totalProgress);
    progress = 0;
    progressBar.show(0);
    // For saving stacks
    int blockSize = 10;
    int nPerThread = (int) Math.ceil((double) frames / nThreads);
    // Convert to fit the block size
    nPerThread = (int) Math.ceil((double) nPerThread / blockSize) * blockSize;
    long seed = start;
    for (int p : photons) {
        statusLine.setText("Simulating " + Utils.pleural(p, "photon"));
        // Create the directory
        File out = new File(directory, String.format("photon%03d", p));
        if (!out.exists())
            out.mkdir();
        for (int from = 0; from < frames; ) {
            int to = Math.min(from + nPerThread, frames);
            futures.add(executor.submit(new SimulationWorker(seed++, out.getPath(), simulationStack, from, to, blockSize, p)));
            from = to;
        }
        // Wait for all to finish
        for (int t = futures.size(); t-- > 0; ) {
            try {
                // The future .get() method will block until completed
                futures.get(t).get();
            } catch (Exception e) {
                // This should not happen. 
                e.printStackTrace();
            }
        }
        futures.clear();
    }
    Utils.setShowStatus(true);
    Utils.setShowProgress(true);
    IJ.showProgress(1);
    executor.shutdown();
    Utils.log("Simulation time = " + Utils.timeToString(System.currentTimeMillis() - start));
}
Also used : PoissonDistribution(org.apache.commons.math3.distribution.PoissonDistribution) TurboList(gdsc.core.utils.TurboList) ImageStack(ij.ImageStack) ExponentialDistribution(org.apache.commons.math3.distribution.ExponentialDistribution) JLabel(javax.swing.JLabel) Well19937c(org.apache.commons.math3.random.Well19937c) ImagePlus(ij.ImagePlus) PseudoRandomGenerator(gdsc.core.utils.PseudoRandomGenerator) RandomGenerator(org.apache.commons.math3.random.RandomGenerator) ExecutorService(java.util.concurrent.ExecutorService) Future(java.util.concurrent.Future) File(java.io.File)

Example 4 with Variance

use of org.apache.commons.math3.stat.descriptive.moment.Variance in project GDSC-SMLM by aherbert.

the class CMOSAnalysis method computeError.

private void computeError(int slice, ImageStack simulationStack) {
    String label = simulationStack.getSliceLabel(slice);
    float[] e = (float[]) simulationStack.getPixels(slice);
    float[] o = (float[]) measuredStack.getPixels(slice);
    // Get the mean error
    Statistics s = new Statistics();
    for (int i = e.length; i-- > 0; ) s.add(o[i] - e[i]);
    StringBuilder result = new StringBuilder("Error ").append(label);
    result.append(" = ").append(Utils.rounded(s.getMean()));
    result.append(" +/- ").append(Utils.rounded(s.getStandardDeviation()));
    // Do statistical tests
    double[] x = Utils.toDouble(e), y = Utils.toDouble(o);
    PearsonsCorrelation c = new PearsonsCorrelation();
    result.append(" : R=").append(Utils.rounded(c.correlation(x, y)));
    // Mann-Whitney U is valid for any distribution, e.g. variance
    MannWhitneyUTest test = new MannWhitneyUTest();
    double p = test.mannWhitneyUTest(x, y);
    result.append(" : Mann-Whitney U p=").append(Utils.rounded(p)).append(' ').append(((p < 0.05) ? "reject" : "accept"));
    if (slice != 2) {
        // T-Test is valid for approximately Normal distributions, e.g. offset and gain
        p = TestUtils.tTest(x, y);
        result.append(" : T-Test p=").append(Utils.rounded(p)).append(' ').append(((p < 0.05) ? "reject" : "accept"));
        p = TestUtils.pairedTTest(x, y);
        result.append(" : Paired T-Test p=").append(Utils.rounded(p)).append(' ').append(((p < 0.05) ? "reject" : "accept"));
    }
    Utils.log(result.toString());
}
Also used : MannWhitneyUTest(org.apache.commons.math3.stat.inference.MannWhitneyUTest) Statistics(gdsc.core.utils.Statistics) PearsonsCorrelation(org.apache.commons.math3.stat.correlation.PearsonsCorrelation)

Example 5 with Variance

use of org.apache.commons.math3.stat.descriptive.moment.Variance in project gatk by broadinstitute.

the class SliceSamplerUnitTest method testSliceSamplingOfPeakedBetaDistribution.

/**
     * Test slice sampling of a peaked beta distribution as an example of sampling of a bounded random variable.
     * Checks that input mean and variance are recovered by 10000 samples to a relative error of 0.5% and 2%,
     * respectively.
     */
@Test
public void testSliceSamplingOfPeakedBetaDistribution() {
    rng.setSeed(RANDOM_SEED);
    final double alpha = 10.;
    final double beta = 4.;
    final BetaDistribution betaDistribution = new BetaDistribution(alpha, beta);
    final Function<Double, Double> betaLogPDF = betaDistribution::logDensity;
    final double xInitial = 0.5;
    final double xMin = 0.;
    final double xMax = 1.;
    final double width = 0.1;
    final int numSamples = 10000;
    final SliceSampler betaSampler = new SliceSampler(rng, betaLogPDF, xMin, xMax, width);
    final double[] samples = Doubles.toArray(betaSampler.sample(xInitial, numSamples));
    final double mean = betaDistribution.getNumericalMean();
    final double variance = betaDistribution.getNumericalVariance();
    final double sampleMean = new Mean().evaluate(samples);
    final double sampleVariance = new Variance().evaluate(samples);
    Assert.assertEquals(relativeError(sampleMean, mean), 0., 0.005);
    Assert.assertEquals(relativeError(sampleVariance, variance), 0., 0.02);
}
Also used : BetaDistribution(org.apache.commons.math3.distribution.BetaDistribution) Mean(org.apache.commons.math3.stat.descriptive.moment.Mean) Variance(org.apache.commons.math3.stat.descriptive.moment.Variance) Test(org.testng.annotations.Test)

Aggregations

Collectors (java.util.stream.Collectors)24 IntStream (java.util.stream.IntStream)24 ParamUtils (org.broadinstitute.hellbender.utils.param.ParamUtils)22 Nonnull (javax.annotation.Nonnull)20 RandomGenerator (org.apache.commons.math3.random.RandomGenerator)18 Variance (org.apache.commons.math3.stat.descriptive.moment.Variance)17 List (java.util.List)16 FastMath (org.apache.commons.math3.util.FastMath)16 Utils (org.broadinstitute.hellbender.utils.Utils)16 INDArray (org.nd4j.linalg.api.ndarray.INDArray)16 Arrays (java.util.Arrays)14 Function (java.util.function.Function)14 Nullable (javax.annotation.Nullable)14 ImmutablePair (org.apache.commons.lang3.tuple.ImmutablePair)14 Logger (org.apache.logging.log4j.Logger)14 GATKException (org.broadinstitute.hellbender.exceptions.GATKException)14 UserException (org.broadinstitute.hellbender.exceptions.UserException)14 Nd4j (org.nd4j.linalg.factory.Nd4j)14 NDArrayIndex (org.nd4j.linalg.indexing.NDArrayIndex)14 VisibleForTesting (com.google.common.annotations.VisibleForTesting)12