Examples with Statistics gdsc.core.utils.Statistics used on opensource projects

Example 1 with Statistics

use of gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.

the class Gaussian2DFunctionTest method functionComputesTargetGradient.

private void functionComputesTargetGradient(int targetParameter) {
    int gradientIndex = findGradientIndex(f1, targetParameter);
    double[] dyda = new double[f1.gradientIndices().length];
    double[] dyda2 = new double[dyda.length];
    double[] a;
    Gaussian2DFunction f1a = GaussianFunctionFactory.create2D(1, maxx, maxy, flags, zModel);
    Gaussian2DFunction f1b = GaussianFunctionFactory.create2D(1, maxx, maxy, flags, zModel);
    Statistics s = new Statistics();
    for (double background : testbackground) // Peak 1
    for (double amplitude1 : testamplitude1) for (double shape1 : testshape1) for (double cx1 : testcx1) for (double cy1 : testcy1) for (double[] w1 : testw1) {
        a = createParameters(background, amplitude1, shape1, cx1, cy1, w1[0], w1[1]);
        f1.initialise(a);
        // Numerically solve gradient. 
        // Calculate the step size h to be an exact numerical representation
        final double xx = a[targetParameter];
        // Get h to minimise roundoff error
        double h = Precision.representableDelta(xx, h_);
        // Evaluate at (x+h) and (x-h)
        a = createParameters(background, amplitude1, shape1, cx1, cy1, w1[0], w1[1]);
        a[targetParameter] = xx + h;
        f1a.initialise(a);
        a = createParameters(background, amplitude1, shape1, cx1, cy1, w1[0], w1[1]);
        a[targetParameter] = xx - h;
        f1b.initialise(a);
        for (int x : testx) for (int y : testy) {
            int i = y * maxx + x;
            f1.eval(i, dyda);
            double value2 = f1a.eval(i, dyda2);
            double value3 = f1b.eval(i, dyda2);
            double gradient = (value2 - value3) / (2 * h);
            double error = DoubleEquality.relativeError(gradient, dyda2[gradientIndex]);
            s.add(error);
            Assert.assertTrue(gradient + " sign != " + dyda2[gradientIndex], (gradient * dyda2[gradientIndex]) >= 0);
            //System.out.printf("[%d,%d] %f == [%d] %f? (%g)\n", x, y, gradient,
            //		gradientIndex, dyda2[gradientIndex], error);
            //System.out.printf("[%d,%d] %f == [%d] %f?\n", x, y, gradient, gradientIndex, dyda[gradientIndex]);
            Assert.assertTrue(gradient + " != " + dyda[gradientIndex], eq.almostEqualRelativeOrAbsolute(gradient, dyda[gradientIndex]));
        }
    }
    System.out.printf("functionComputesTargetGradient %s %s (error %s +/- %s)\n", f1.getClass().getSimpleName(), f1.getName(targetParameter), Utils.rounded(s.getMean()), Utils.rounded(s.getStandardDeviation()));
}

Example 2 with Statistics

use of gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.

the class DiffusionRateTest method run.

/*
	 * (non-Javadoc)
	 * 
	 * @see ij.plugin.PlugIn#run(java.lang.String)
	 */
public void run(String arg) {
    SMLMUsageTracker.recordPlugin(this.getClass(), arg);
    if (IJ.controlKeyDown()) {
        simpleTest();
        return;
    }
    extraOptions = Utils.isExtraOptions();
    if (!showDialog())
        return;
    lastSimulatedDataset[0] = lastSimulatedDataset[1] = "";
    lastSimulatedPrecision = 0;
    final int totalSteps = (int) Math.ceil(settings.seconds * settings.stepsPerSecond);
    conversionFactor = 1000000.0 / (settings.pixelPitch * settings.pixelPitch);
    // Diffusion rate is um^2/sec. Convert to pixels per simulation frame.
    final double diffusionRateInPixelsPerSecond = settings.diffusionRate * conversionFactor;
    final double diffusionRateInPixelsPerStep = diffusionRateInPixelsPerSecond / settings.stepsPerSecond;
    final double precisionInPixels = myPrecision / settings.pixelPitch;
    final boolean addError = myPrecision != 0;
    Utils.log(TITLE + " : D = %s um^2/sec, Precision = %s nm", Utils.rounded(settings.diffusionRate, 4), Utils.rounded(myPrecision, 4));
    Utils.log("Mean-displacement per dimension = %s nm/sec", Utils.rounded(1e3 * ImageModel.getRandomMoveDistance(settings.diffusionRate), 4));
    if (extraOptions)
        Utils.log("Step size = %s, precision = %s", Utils.rounded(ImageModel.getRandomMoveDistance(diffusionRateInPixelsPerStep)), Utils.rounded(precisionInPixels));
    // Convert diffusion co-efficient into the standard deviation for the random walk
    final double diffusionSigma = (settings.getDiffusionType() == DiffusionType.LINEAR_WALK) ? // Q. What should this be? At the moment just do 1D diffusion on a random vector
    ImageModel.getRandomMoveDistance(diffusionRateInPixelsPerStep) : ImageModel.getRandomMoveDistance(diffusionRateInPixelsPerStep);
    Utils.log("Simulation step-size = %s nm", Utils.rounded(settings.pixelPitch * diffusionSigma, 4));
    // Move the molecules and get the diffusion rate
    IJ.showStatus("Simulating ...");
    final long start = System.nanoTime();
    final long seed = System.currentTimeMillis() + System.identityHashCode(this);
    RandomGenerator[] random = new RandomGenerator[3];
    RandomGenerator[] random2 = new RandomGenerator[3];
    for (int i = 0; i < 3; i++) {
        random[i] = new Well19937c(seed + i * 12436);
        random2[i] = new Well19937c(seed + i * 678678 + 3);
    }
    Statistics[] stats2D = new Statistics[totalSteps];
    Statistics[] stats3D = new Statistics[totalSteps];
    StoredDataStatistics jumpDistances2D = new StoredDataStatistics(totalSteps);
    StoredDataStatistics jumpDistances3D = new StoredDataStatistics(totalSteps);
    for (int j = 0; j < totalSteps; j++) {
        stats2D[j] = new Statistics();
        stats3D[j] = new Statistics();
    }
    SphericalDistribution dist = new SphericalDistribution(settings.confinementRadius / settings.pixelPitch);
    Statistics asymptote = new Statistics();
    // Save results to memory
    MemoryPeakResults results = new MemoryPeakResults(totalSteps);
    Calibration cal = new Calibration(settings.pixelPitch, 1, 1000.0 / settings.stepsPerSecond);
    results.setCalibration(cal);
    results.setName(TITLE);
    int peak = 0;
    // Store raw coordinates
    ArrayList<Point> points = new ArrayList<Point>(totalSteps);
    StoredData totalJumpDistances1D = new StoredData(settings.particles);
    StoredData totalJumpDistances2D = new StoredData(settings.particles);
    StoredData totalJumpDistances3D = new StoredData(settings.particles);
    for (int i = 0; i < settings.particles; i++) {
        if (i % 16 == 0) {
            IJ.showProgress(i, settings.particles);
            if (Utils.isInterrupted())
                return;
        }
        // Increment the frame so that tracing analysis can distinguish traces
        peak++;
        double[] origin = new double[3];
        final int id = i + 1;
        MoleculeModel m = new MoleculeModel(id, origin.clone());
        if (addError)
            origin = addError(origin, precisionInPixels, random);
        if (useConfinement) {
            // Note: When using confinement the average displacement should asymptote
            // at the average distance of a point from the centre of a ball. This is 3r/4.
            // See: http://answers.yahoo.com/question/index?qid=20090131162630AAMTUfM
            // The equivalent in 2D is 2r/3. However although we are plotting 2D distance
            // this is a projection of the 3D position onto the plane and so the particles
            // will not be evenly spread (there will be clustering at centre caused by the
            // poles)
            final double[] axis = (settings.getDiffusionType() == DiffusionType.LINEAR_WALK) ? nextVector() : null;
            for (int j = 0; j < totalSteps; j++) {
                double[] xyz = m.getCoordinates();
                double[] originalXyz = xyz.clone();
                for (int n = confinementAttempts; n-- > 0; ) {
                    if (settings.getDiffusionType() == DiffusionType.GRID_WALK)
                        m.walk(diffusionSigma, random);
                    else if (settings.getDiffusionType() == DiffusionType.LINEAR_WALK)
                        m.slide(diffusionSigma, axis, random[0]);
                    else
                        m.move(diffusionSigma, random);
                    if (!dist.isWithin(m.getCoordinates())) {
                        // Reset position
                        for (int k = 0; k < 3; k++) xyz[k] = originalXyz[k];
                    } else {
                        // The move was allowed
                        break;
                    }
                }
                points.add(new Point(id, xyz));
                if (addError)
                    xyz = addError(xyz, precisionInPixels, random2);
                peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
            }
            asymptote.add(distance(m.getCoordinates()));
        } else {
            if (settings.getDiffusionType() == DiffusionType.GRID_WALK) {
                for (int j = 0; j < totalSteps; j++) {
                    m.walk(diffusionSigma, random);
                    double[] xyz = m.getCoordinates();
                    points.add(new Point(id, xyz));
                    if (addError)
                        xyz = addError(xyz, precisionInPixels, random2);
                    peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
                }
            } else if (settings.getDiffusionType() == DiffusionType.LINEAR_WALK) {
                final double[] axis = nextVector();
                for (int j = 0; j < totalSteps; j++) {
                    m.slide(diffusionSigma, axis, random[0]);
                    double[] xyz = m.getCoordinates();
                    points.add(new Point(id, xyz));
                    if (addError)
                        xyz = addError(xyz, precisionInPixels, random2);
                    peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
                }
            } else {
                for (int j = 0; j < totalSteps; j++) {
                    m.move(diffusionSigma, random);
                    double[] xyz = m.getCoordinates();
                    points.add(new Point(id, xyz));
                    if (addError)
                        xyz = addError(xyz, precisionInPixels, random2);
                    peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
                }
            }
        }
        // Debug: record all the particles so they can be analysed
        // System.out.printf("%f %f %f\n", m.getX(), m.getY(), m.getZ());
        final double[] xyz = m.getCoordinates();
        double d2 = 0;
        totalJumpDistances1D.add(d2 = xyz[0] * xyz[0]);
        totalJumpDistances2D.add(d2 += xyz[1] * xyz[1]);
        totalJumpDistances3D.add(d2 += xyz[2] * xyz[2]);
    }
    final double time = (System.nanoTime() - start) / 1000000.0;
    IJ.showProgress(1);
    MemoryPeakResults.addResults(results);
    lastSimulatedDataset[0] = results.getName();
    lastSimulatedPrecision = myPrecision;
    // Convert pixels^2/step to um^2/sec
    final double msd2D = (jumpDistances2D.getMean() / conversionFactor) / (results.getCalibration().getExposureTime() / 1000);
    final double msd3D = (jumpDistances3D.getMean() / conversionFactor) / (results.getCalibration().getExposureTime() / 1000);
    Utils.log("Raw data D=%s um^2/s, Precision = %s nm, N=%d, step=%s s, mean2D=%s um^2, MSD 2D = %s um^2/s, mean3D=%s um^2, MSD 3D = %s um^2/s", Utils.rounded(settings.diffusionRate), Utils.rounded(myPrecision), jumpDistances2D.getN(), Utils.rounded(results.getCalibration().getExposureTime() / 1000), Utils.rounded(jumpDistances2D.getMean() / conversionFactor), Utils.rounded(msd2D), Utils.rounded(jumpDistances3D.getMean() / conversionFactor), Utils.rounded(msd3D));
    aggregateIntoFrames(points, addError, precisionInPixels, random2);
    IJ.showStatus("Analysing results ...");
    if (showDiffusionExample) {
        showExample(totalSteps, diffusionSigma, random);
    }
    // Plot a graph of mean squared distance
    double[] xValues = new double[stats2D.length];
    double[] yValues2D = new double[stats2D.length];
    double[] yValues3D = new double[stats3D.length];
    double[] upper2D = new double[stats2D.length];
    double[] lower2D = new double[stats2D.length];
    double[] upper3D = new double[stats3D.length];
    double[] lower3D = new double[stats3D.length];
    SimpleRegression r2D = new SimpleRegression(false);
    SimpleRegression r3D = new SimpleRegression(false);
    final int firstN = (useConfinement) ? fitN : totalSteps;
    for (int j = 0; j < totalSteps; j++) {
        // Convert steps to seconds
        xValues[j] = (double) (j + 1) / settings.stepsPerSecond;
        // Convert values in pixels^2 to um^2
        final double mean2D = stats2D[j].getMean() / conversionFactor;
        final double mean3D = stats3D[j].getMean() / conversionFactor;
        final double sd2D = stats2D[j].getStandardDeviation() / conversionFactor;
        final double sd3D = stats3D[j].getStandardDeviation() / conversionFactor;
        yValues2D[j] = mean2D;
        yValues3D[j] = mean3D;
        upper2D[j] = mean2D + sd2D;
        lower2D[j] = mean2D - sd2D;
        upper3D[j] = mean3D + sd3D;
        lower3D[j] = mean3D - sd3D;
        if (j < firstN) {
            r2D.addData(xValues[j], yValues2D[j]);
            r3D.addData(xValues[j], yValues3D[j]);
        }
    }
    // TODO - Fit using the equation for 2D confined diffusion:
    // MSD = 4s^2 + R^2 (1 - 0.99e^(-1.84^2 Dt / R^2)
    // s = localisation precision
    // R = confinement radius
    // D = 2D diffusion coefficient
    // t = time
    final PolynomialFunction fitted2D, fitted3D;
    if (r2D.getN() > 0) {
        // Do linear regression to get diffusion rate
        final double[] best2D = new double[] { r2D.getIntercept(), r2D.getSlope() };
        fitted2D = new PolynomialFunction(best2D);
        final double[] best3D = new double[] { r3D.getIntercept(), r3D.getSlope() };
        fitted3D = new PolynomialFunction(best3D);
        // For 2D diffusion: d^2 = 4D
        // where: d^2 = mean-square displacement
        double D = best2D[1] / 4.0;
        String msg = "2D Diffusion rate = " + Utils.rounded(D, 4) + " um^2 / sec (" + Utils.timeToString(time) + ")";
        IJ.showStatus(msg);
        Utils.log(msg);
        D = best3D[1] / 6.0;
        Utils.log("3D Diffusion rate = " + Utils.rounded(D, 4) + " um^2 / sec (" + Utils.timeToString(time) + ")");
    } else {
        fitted2D = fitted3D = null;
    }
    // Create plots
    plotMSD(totalSteps, xValues, yValues2D, lower2D, upper2D, fitted2D, 2);
    plotMSD(totalSteps, xValues, yValues3D, lower3D, upper3D, fitted3D, 3);
    plotJumpDistances(TITLE, jumpDistances2D, 2, 1);
    plotJumpDistances(TITLE, jumpDistances3D, 3, 1);
    if (idCount > 0)
        new WindowOrganiser().tileWindows(idList);
    if (useConfinement)
        Utils.log("3D asymptote distance = %s nm (expected %.2f)", Utils.rounded(asymptote.getMean() * settings.pixelPitch, 4), 3 * settings.confinementRadius / 4);
}

Example 3 with Statistics

use of gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.

the class DataEstimator method getEstimate.

private void getEstimate() {
    if (estimate == null) {
        estimate = new float[5];
        if (h == null) {
            h = FloatHistogram.buildHistogram(data.clone(), true);
            h = h.compact(histogramSize);
        }
        // Threshold the data
        final float t = estimate[ESTIMATE_THRESHOLD] = h.getAutoThreshold(thresholdMethod);
        // Get stats below the threshold
        Statistics stats = new Statistics();
        for (int i = h.minBin; i <= h.maxBin; i++) {
            if (h.getValue(i) >= t)
                break;
            stats.add(h.h[i], h.getValue(i));
        }
        // Check if background region is large enough
        estimate[ESTIMATE_BACKGROUND_SIZE] = stats.getN();
        if (stats.getN() > fraction * data.length) {
            // Background region is large enough
            estimate[ESTIMATE_LARGE_ENOUGH] = 1;
        } else {
            // Recompute with all the data
            stats = new Statistics(data);
        }
        estimate[ESTIMATE_BACKGROUND] = (float) stats.getMean();
        estimate[ESTIMATE_NOISE] = (float) stats.getStandardDeviation();
    }
}

Example 4 with Statistics

use of gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.

the class TraceDiffusion method run.

/*
	 * (non-Javadoc)
	 * 
	 * @see ij.plugin.PlugIn#run(java.lang.String)
	 */
public void run(String arg) {
    SMLMUsageTracker.recordPlugin(this.getClass(), arg);
    jumpDistanceParameters = null;
    extraOptions = Utils.isExtraOptions();
    if (MemoryPeakResults.isMemoryEmpty()) {
        IJ.error(TITLE, "No localisations in memory");
        return;
    }
    ArrayList<MemoryPeakResults> allResults = new ArrayList<MemoryPeakResults>();
    // Option to pick multiple input datasets together using a list box.
    if ("multi".equals(arg)) {
        if (!showMultiDialog(allResults))
            return;
    }
    // This shows the dialog for selecting trace options
    if (!showTraceDialog(allResults))
        return;
    if (// Sense check
    allResults.isEmpty())
        return;
    Utils.log(TITLE + "...");
    // This optionally collects additional datasets then gets the traces:
    // - Trace each single dataset (and store in memory)
    // - Combine trace results held in memory
    Trace[] traces = getTraces(allResults);
    // This still allows a zero entry in the results table.
    if (traces.length > 0)
        if (!showDialog())
            return;
    int count = traces.length;
    double[] fitMSDResult = null;
    int n = 0;
    double[][] jdParams = null;
    if (count > 0) {
        calculatePrecision(traces, allResults.size() > 1);
        //--- MSD Analysis ---
        // Conversion constants
        final double px2ToUm2 = results.getCalibration().getNmPerPixel() * results.getCalibration().getNmPerPixel() / 1e6;
        final double px2ToUm2PerSecond = px2ToUm2 / exposureTime;
        // Get the maximum trace length
        int length = settings.minimumTraceLength;
        if (!settings.truncate) {
            for (Trace trace : traces) {
                if (length < trace.size())
                    length = trace.size();
            }
        }
        // Get the localisation error (4s^2) in um^2
        final double error = (settings.precisionCorrection) ? 4 * precision * precision / 1e6 : 0;
        // Pre-calculate MSD correction factors. This accounts for the fact that the distance moved 
        // in the start/end frames is reduced due to the averaging of the particle location over the 
        // entire frame into a single point. The true MSD may be restored by applying a factor.
        // Note: These are used for the calculation of the diffusion coefficients per molecule and 
        // the MSD passed to the Jump Distance analysis. However the error is not included in the 
        // jump distance analysis so will be subtracted from the fitted D coefficients later.
        final double[] factors;
        if (settings.msdCorrection) {
            factors = new double[length];
            for (int t = 1; t < length; t++) factors[t] = JumpDistanceAnalysis.getConversionfactor(t);
        } else {
            factors = Utils.newArray(length, 0.0, 1.0);
        }
        // Extract the mean-squared distance statistics
        Statistics[] stats = new Statistics[length];
        for (int i = 0; i < stats.length; i++) stats[i] = new Statistics();
        ArrayList<double[]> distances = (saveTraceDistances || displayTraceLength) ? new ArrayList<double[]>(traces.length) : null;
        // Store all the jump distances at the specified interval
        StoredDataStatistics jumpDistances = new StoredDataStatistics();
        final int jumpDistanceInterval = settings.jumpDistance;
        // Compute squared distances
        StoredDataStatistics msdPerMoleculeAllVsAll = new StoredDataStatistics();
        StoredDataStatistics msdPerMoleculeAdjacent = new StoredDataStatistics();
        for (Trace trace : traces) {
            ArrayList<PeakResult> results = trace.getPoints();
            // Sum the MSD and the time
            final int traceLength = (settings.truncate) ? settings.minimumTraceLength : trace.size();
            // Get the mean for each time separation
            double[] sumDistance = new double[traceLength + 1];
            double[] sumTime = new double[sumDistance.length];
            // Do the distances to the origin (saving if necessary)
            {
                final float x = results.get(0).getXPosition();
                final float y = results.get(0).getYPosition();
                if (distances != null) {
                    double[] msd = new double[traceLength - 1];
                    for (int j = 1; j < traceLength; j++) {
                        final int t = j;
                        final double d = distance2(x, y, results.get(j));
                        msd[j - 1] = px2ToUm2 * d;
                        if (t == jumpDistanceInterval)
                            jumpDistances.add(msd[j - 1]);
                        sumDistance[t] += d;
                        sumTime[t] += t;
                    }
                    distances.add(msd);
                } else {
                    for (int j = 1; j < traceLength; j++) {
                        final int t = j;
                        final double d = distance2(x, y, results.get(j));
                        if (t == jumpDistanceInterval)
                            jumpDistances.add(px2ToUm2 * d);
                        sumDistance[t] += d;
                        sumTime[t] += t;
                    }
                }
            }
            if (settings.internalDistances) {
                // Do the internal distances
                for (int i = 1; i < traceLength; i++) {
                    final float x = results.get(i).getXPosition();
                    final float y = results.get(i).getYPosition();
                    for (int j = i + 1; j < traceLength; j++) {
                        final int t = j - i;
                        final double d = distance2(x, y, results.get(j));
                        if (t == jumpDistanceInterval)
                            jumpDistances.add(px2ToUm2 * d);
                        sumDistance[t] += d;
                        sumTime[t] += t;
                    }
                }
                // Add the average distance per time separation to the population
                for (int t = 1; t < traceLength; t++) {
                    // Note: (traceLength - t) == count
                    stats[t].add(sumDistance[t] / (traceLength - t));
                }
            } else {
                // Add the distance per time separation to the population
                for (int t = 1; t < traceLength; t++) {
                    stats[t].add(sumDistance[t]);
                }
            }
            // Fix this for the precision and MSD adjustment.
            // It may be necessary to:
            // - sum the raw distances for each time interval (this is sumDistance[t])
            // - subtract the precision error
            // - apply correction factor for the n-frames to get actual MSD
            // - sum the actual MSD
            double sumD = 0, sumD_adjacent = Math.max(0, sumDistance[1] - error) * factors[1];
            double sumT = 0, sumT_adjacent = sumTime[1];
            for (int t = 1; t < traceLength; t++) {
                sumD += Math.max(0, sumDistance[t] - error) * factors[t];
                sumT += sumTime[t];
            }
            // Calculate the average displacement for the trace (do not simply use the largest 
            // time separation since this will miss moving molecules that end up at the origin)
            msdPerMoleculeAllVsAll.add(px2ToUm2PerSecond * sumD / sumT);
            msdPerMoleculeAdjacent.add(px2ToUm2PerSecond * sumD_adjacent / sumT_adjacent);
        }
        StoredDataStatistics dPerMoleculeAllVsAll = null;
        StoredDataStatistics dPerMoleculeAdjacent = null;
        if (saveTraceDistances || (settings.showHistograms && displayDHistogram)) {
            dPerMoleculeAllVsAll = calculateDiffusionCoefficient(msdPerMoleculeAllVsAll);
            dPerMoleculeAdjacent = calculateDiffusionCoefficient(msdPerMoleculeAdjacent);
        }
        if (saveTraceDistances) {
            saveTraceDistances(traces.length, distances, msdPerMoleculeAllVsAll, msdPerMoleculeAdjacent, dPerMoleculeAllVsAll, dPerMoleculeAdjacent);
        }
        if (displayTraceLength) {
            StoredDataStatistics lengths = calculateTraceLengths(distances);
            showHistogram(lengths, "Trace length (um)");
        }
        if (displayTraceSize) {
            StoredDataStatistics sizes = calculateTraceSizes(traces);
            showHistogram(sizes, "Trace size", true);
        }
        // Plot the per-trace histogram of MSD and D
        if (settings.showHistograms) {
            if (displayMSDHistogram) {
                showHistogram(msdPerMoleculeAllVsAll, "MSD/Molecule (all-vs-all)");
                showHistogram(msdPerMoleculeAdjacent, "MSD/Molecule (adjacent)");
            }
            if (displayDHistogram) {
                showHistogram(dPerMoleculeAllVsAll, "D/Molecule (all-vs-all)");
                showHistogram(dPerMoleculeAdjacent, "D/Molecule (adjacent)");
            }
        }
        // Calculate the mean squared distance (MSD)
        double[] x = new double[stats.length];
        double[] y = new double[x.length];
        double[] sd = new double[x.length];
        // Intercept is the 4s^2 (in um^2)
        y[0] = 4 * precision * precision / 1e6;
        for (int i = 1; i < stats.length; i++) {
            x[i] = i * exposureTime;
            y[i] = stats[i].getMean() * px2ToUm2;
            //sd[i] = stats[i].getStandardDeviation() * px2ToUm2;
            sd[i] = stats[i].getStandardError() * px2ToUm2;
        }
        String title = TITLE + " MSD";
        Plot2 plot = plotMSD(x, y, sd, title);
        // Fit the MSD using a linear fit
        fitMSDResult = fitMSD(x, y, title, plot);
        // Jump Distance analysis
        if (saveRawData)
            saveStatistics(jumpDistances, "Jump Distance", "Distance (um^2)", false);
        // Calculate the cumulative jump-distance histogram
        double[][] jdHistogram = JumpDistanceAnalysis.cumulativeHistogram(jumpDistances.getValues());
        // Always show the jump distance histogram
        jdTitle = TITLE + " Jump Distance";
        jdPlot = new Plot2(jdTitle, "Distance (um^2)", "Cumulative Probability", jdHistogram[0], jdHistogram[1]);
        display(jdTitle, jdPlot);
        // Fit Jump Distance cumulative probability
        n = jumpDistances.getN();
        jumpDistanceParameters = jdParams = fitJumpDistance(jumpDistances, jdHistogram);
    }
    summarise(traces, fitMSDResult, n, jdParams);
}

Example 5 with Statistics

use of gdsc.core.utils.Statistics 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());
}