use of org.apache.commons.math3.stat.descriptive.rank.Min in project ffx by mjschnie.
the class ReflectionList method ordinal.
/**
* <p>
* ordinal</p>
*
* @param s a double.
* @return a double.
*/
public double ordinal(double s) {
double r = (s - minres) / (maxres - minres);
r = min(r, 0.999) * 1000.0;
int i = (int) r;
r -= floor(r);
return ((1.0 - r) * hist[i] + r * hist[i + 1]);
}
use of org.apache.commons.math3.stat.descriptive.rank.Min in project GDSC-SMLM by aherbert.
the class CreateData method createUniformDistribution.
/**
* Create distribution within an XY border
*
* @param border
* @return
*/
private UniformDistribution createUniformDistribution(double border) {
double depth = (settings.fixedDepth) ? settings.depth / settings.pixelPitch : settings.depth / (2 * settings.pixelPitch);
// Ensure the focal plane is in the middle of the zDepth
double[] max = new double[] { settings.size / 2 - border, settings.size / 2 - border, depth };
double[] min = new double[3];
for (int i = 0; i < 3; i++) min[i] = -max[i];
if (settings.fixedDepth)
min[2] = max[2];
// Try using different distributions:
final RandomGenerator rand1 = createRandomGenerator();
if (settings.distribution.equals(DISTRIBUTION[UNIFORM_HALTON])) {
return new UniformDistribution(min, max, rand1.nextInt());
}
if (settings.distribution.equals(DISTRIBUTION[UNIFORM_SOBOL])) {
SobolSequenceGenerator rvg = new SobolSequenceGenerator(3);
rvg.skipTo(rand1.nextInt());
return new UniformDistribution(min, max, rvg);
}
// Create a distribution using random generators for each dimension
UniformDistribution distribution = new UniformDistribution(min, max, this);
return distribution;
}
use of org.apache.commons.math3.stat.descriptive.rank.Min in project GDSC-SMLM by aherbert.
the class PCPALMMolecules method calculateAveragePrecision.
/**
* Calculate the average precision by fitting a skewed Gaussian to the histogram of the precision distribution.
* <p>
* A simple mean and SD of the histogram is computed. If the mean of the Skewed Gaussian does not fit within 3 SDs
* of the simple mean then the simple mean is returned.
*
* @param molecules
* @param title
* the plot title (null if no plot should be displayed)
* @param histogramBins
* @param logFitParameters
* Record the fit parameters to the ImageJ log
* @param removeOutliers
* The distribution is created using all values within 1.5x the inter-quartile range (IQR) of the data
* @return The average precision
*/
public double calculateAveragePrecision(ArrayList<Molecule> molecules, String title, int histogramBins, boolean logFitParameters, boolean removeOutliers) {
// Plot histogram of the precision
float[] data = new float[molecules.size()];
DescriptiveStatistics stats = new DescriptiveStatistics();
double yMin = Double.NEGATIVE_INFINITY, yMax = 0;
for (int i = 0; i < data.length; i++) {
data[i] = (float) molecules.get(i).precision;
stats.addValue(data[i]);
}
// Set the min and max y-values using 1.5 x IQR
if (removeOutliers) {
double lower = stats.getPercentile(25);
double upper = stats.getPercentile(75);
if (Double.isNaN(lower) || Double.isNaN(upper)) {
if (logFitParameters)
Utils.log("Error computing IQR: %f - %f", lower, upper);
} else {
double iqr = upper - lower;
yMin = FastMath.max(lower - iqr, stats.getMin());
yMax = FastMath.min(upper + iqr, stats.getMax());
if (logFitParameters)
Utils.log(" Data range: %f - %f. Plotting 1.5x IQR: %f - %f", stats.getMin(), stats.getMax(), yMin, yMax);
}
}
if (yMin == Double.NEGATIVE_INFINITY) {
yMin = stats.getMin();
yMax = stats.getMax();
if (logFitParameters)
Utils.log(" Data range: %f - %f", yMin, yMax);
}
float[][] hist = Utils.calcHistogram(data, yMin, yMax, histogramBins);
Plot2 plot = null;
if (title != null) {
plot = new Plot2(title, "Precision", "Frequency");
float[] xValues = hist[0];
float[] yValues = hist[1];
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, Maths.max(yValues) * 1.05);
}
plot.addPoints(xValues, yValues, Plot2.BAR);
Utils.display(title, plot);
}
// Extract non-zero data
float[] x = Arrays.copyOf(hist[0], hist[0].length);
float[] y = hist[1];
int count = 0;
float dx = (x[1] - x[0]) * 0.5f;
for (int i = 0; i < y.length; i++) if (y[i] > 0) {
x[count] = x[i] + dx;
y[count] = y[i];
count++;
}
x = Arrays.copyOf(x, count);
y = Arrays.copyOf(y, count);
// Sense check to fitted data. Get mean and SD of histogram
double[] stats2 = Utils.getHistogramStatistics(x, y);
double mean = stats2[0];
if (logFitParameters)
log(" Initial Statistics: %f +/- %f", stats2[0], stats2[1]);
// Standard Gaussian fit
double[] parameters = fitGaussian(x, y);
if (parameters == null) {
log(" Failed to fit initial Gaussian");
return mean;
}
double newMean = parameters[1];
double error = Math.abs(stats2[0] - newMean) / stats2[1];
if (error > 3) {
log(" Failed to fit Gaussian: %f standard deviations from histogram mean", error);
return mean;
}
if (newMean < yMin || newMean > yMax) {
log(" Failed to fit Gaussian: %f outside data range %f - %f", newMean, yMin, yMax);
return mean;
}
mean = newMean;
if (logFitParameters)
log(" Initial Gaussian: %f @ %f +/- %f", parameters[0], parameters[1], parameters[2]);
double[] initialSolution = new double[] { parameters[0], parameters[1], parameters[2], -1 };
// Fit to a skewed Gaussian (or appropriate function)
double[] skewParameters = fitSkewGaussian(x, y, initialSolution);
if (skewParameters == null) {
log(" Failed to fit Skewed Gaussian");
return mean;
}
SkewNormalFunction sn = new SkewNormalFunction(skewParameters);
if (logFitParameters)
log(" Skewed Gaussian: %f @ %f +/- %f (a = %f) => %f +/- %f", skewParameters[0], skewParameters[1], skewParameters[2], skewParameters[3], sn.getMean(), Math.sqrt(sn.getVariance()));
newMean = sn.getMean();
error = Math.abs(stats2[0] - newMean) / stats2[1];
if (error > 3) {
log(" Failed to fit Skewed Gaussian: %f standard deviations from histogram mean", error);
return mean;
}
if (newMean < yMin || newMean > yMax) {
log(" Failed to fit Skewed Gaussian: %f outside data range %f - %f", newMean, yMin, yMax);
return mean;
}
// Use original histogram x-axis to maintain all the bins
if (plot != null) {
x = hist[0];
for (int i = 0; i < y.length; i++) x[i] += dx;
plot.setColor(Color.red);
addToPlot(plot, x, skewParameters, Plot2.LINE);
plot.setColor(Color.black);
Utils.display(title, plot);
}
// Return the average precision from the fitted curve
return newMean;
}
use of org.apache.commons.math3.stat.descriptive.rank.Min in project GDSC-SMLM by aherbert.
the class PoissonFunctionTest method cumulativeProbabilityIsOneWithRealAbove4.
private void cumulativeProbabilityIsOneWithRealAbove4(final double gain, final double mu, int min, int max) {
final double o = mu * gain;
final PoissonFunction f = new PoissonFunction(1.0 / gain, true);
double p = 0;
SimpsonIntegrator in = new SimpsonIntegrator(3, 30);
try {
p = in.integrate(Integer.MAX_VALUE, new UnivariateFunction() {
public double value(double x) {
return f.likelihood(x, o);
}
}, min, max);
System.out.printf("g=%f, mu=%f, o=%f, p=%f\n", gain, mu, o, p);
Assert.assertEquals(String.format("g=%f, mu=%f", gain, mu), 1, p, 0.02);
} catch (TooManyEvaluationsException e) {
//double inc = max / 20000.0;
//for (double x = 0; x <= max; x += inc)
//{
// final double pp = f.likelihood(x, o);
// //System.out.printf("g=%f, mu=%f, o=%f, p=%f\n", gain, mu, o, pp);
// p += pp;
//}
//System.out.printf("g=%f, mu=%f, o=%f, p=%f\n", gain, mu, o, p);
Assert.assertFalse(e.getMessage(), true);
}
}
use of org.apache.commons.math3.stat.descriptive.rank.Min in project GDSC-SMLM by aherbert.
the class PoissonCalculatorTest method cumulativeProbabilityIsOneWithRealData.
private void cumulativeProbabilityIsOneWithRealData(final double mu, double min, double max, boolean test) {
double p = 0;
UnivariateIntegrator in = new SimpsonIntegrator();
p = in.integrate(20000, new UnivariateFunction() {
public double value(double x) {
double v;
v = PoissonCalculator.likelihood(mu, x);
//System.out.printf("x=%f, v=%f\n", x, v);
return v;
}
}, min, max);
System.out.printf("mu=%f, p=%f\n", mu, p);
if (test) {
Assert.assertEquals(String.format("mu=%f", mu), P_LIMIT, p, 0.02);
}
}
Aggregations