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);
}
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);
}
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));
}
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());
}
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);
}
Aggregations