Examples with FisherInformationMatrix gdsc.smlm.fitting.FisherInformationMatrix used on opensource projects

Example 6 with FisherInformationMatrix

use of gdsc.smlm.fitting.FisherInformationMatrix in project GDSC-SMLM by aherbert.

the class SteppingFunctionSolver method computeDeviations.

/**
	 * Compute the deviations for the parameters a from the last call to
	 * {@link #computeFitValue(double[], double[])}.
	 *
	 * @param a_dev
	 *            the parameter deviations
	 */
protected void computeDeviations(double[] a_dev) {
    // Use a dedicated solver optimised for inverting the matrix diagonal. 
    // The last Hessian matrix should be stored in the working alpha.
    final FisherInformationMatrix m = computeFisherInformationMatrix();
    // This may fail if the matrix cannot be inverted
    final double[] crlb = m.crlb();
    if (crlb == null)
        throw new FunctionSolverException(FitStatus.SINGULAR_NON_LINEAR_SOLUTION);
    setDeviations(a_dev, crlb);
// Use this method for robustness, i.e. it will not fail
//setDeviations(a_dev, m.crlb(true));
}

Example 7 with FisherInformationMatrix

use of gdsc.smlm.fitting.FisherInformationMatrix in project GDSC-SMLM by aherbert.

the class PoissonGradientProcedureTest method crlbIsHigherWithPrecomputed.

@Test
public void crlbIsHigherWithPrecomputed() {
    int iter = 10;
    rdg = new RandomDataGenerator(new Well19937c(30051977));
    ErfGaussian2DFunction func = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(1, 10, 10, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    double[] a = new double[7];
    int n = func.getNumberOfGradients();
    // Get a background
    double[] b = new double[func.size()];
    for (int i = 0; i < b.length; i++) b[i] = rdg.nextUniform(1, 2);
    for (int i = 0; i < iter; i++) {
        a[Gaussian2DFunction.BACKGROUND] = rdg.nextUniform(0.1, 0.3);
        a[Gaussian2DFunction.SIGNAL] = rdg.nextUniform(100, 300);
        a[Gaussian2DFunction.X_POSITION] = rdg.nextUniform(4, 6);
        a[Gaussian2DFunction.Y_POSITION] = rdg.nextUniform(4, 6);
        a[Gaussian2DFunction.X_SD] = rdg.nextUniform(1, 1.3);
        a[Gaussian2DFunction.Y_SD] = rdg.nextUniform(1, 1.3);
        PoissonGradientProcedure p1 = PoissonGradientProcedureFactory.create(func);
        p1.computeFisherInformation(a);
        PoissonGradientProcedure p2 = PoissonGradientProcedureFactory.create(PrecomputedGradient1Function.wrapGradient1Function(func, b));
        p2.computeFisherInformation(a);
        FisherInformationMatrix m1 = new FisherInformationMatrix(p1.getLinear(), n);
        FisherInformationMatrix m2 = new FisherInformationMatrix(p2.getLinear(), n);
        double[] crlb1 = m1.crlb();
        double[] crlb2 = m2.crlb();
        Assert.assertNotNull(crlb1);
        Assert.assertNotNull(crlb2);
        //System.out.printf("%s : %s\n", Arrays.toString(crlb1), Arrays.toString(crlb2));
        for (int j = 0; j < n; j++) Assert.assertTrue(crlb1[j] < crlb2[j]);
    }
}

Example 8 with FisherInformationMatrix

use of gdsc.smlm.fitting.FisherInformationMatrix in project GDSC-SMLM by aherbert.

the class BaseFunctionSolverTest method canFitSingleGaussian.

void canFitSingleGaussian(FunctionSolver solver, boolean applyBounds, NoiseModel noiseModel) {
    // Allow reporting the fit deviations
    boolean report = false;
    double[] crlb = null;
    SimpleArrayMoment m = null;
    double[] noise = getNoise(noiseModel);
    if (solver.isWeighted())
        solver.setWeights(getWeights(noiseModel));
    randomGenerator.setSeed(seed);
    for (double s : signal) {
        double[] expected = createParams(1, s, 0, 0, 1);
        double[] lower = createParams(0, s * 0.5, -0.2, -0.2, 0.8);
        double[] upper = createParams(3, s * 2, 0.2, 0.2, 1.2);
        if (applyBounds)
            solver.setBounds(lower, upper);
        if (report) {
            // Compute the CRLB for a Poisson process
            PoissonGradientProcedure gp = PoissonGradientProcedureFactory.create((Gradient1Function) ((BaseFunctionSolver) solver).getGradientFunction());
            gp.computeFisherInformation(expected);
            FisherInformationMatrix f = new FisherInformationMatrix(gp.getLinear(), gp.n);
            crlb = f.crlbSqrt();
            // Compute the deviations
            m = new SimpleArrayMoment();
        }
        double[] data = drawGaussian(expected, noise, noiseModel);
        for (double db : base) for (double dx : shift) for (double dy : shift) for (double dsx : factor) {
            double[] p = createParams(db, s, dx, dy, dsx);
            double[] fp = fitGaussian(solver, data, p, expected);
            for (int i = 0; i < expected.length; i++) {
                if (fp[i] < lower[i])
                    Assert.assertTrue(String.format("Fit Failed: [%d] %.2f < %.2f: %s != %s", i, fp[i], lower[i], Arrays.toString(fp), Arrays.toString(expected)), false);
                if (fp[i] > upper[i])
                    Assert.assertTrue(String.format("Fit Failed: [%d] %.2f > %.2f: %s != %s", i, fp[i], upper[i], Arrays.toString(fp), Arrays.toString(expected)), false);
                if (report)
                    fp[i] = expected[i] - fp[i];
            }
            // Store the deviations
            if (report)
                m.add(fp);
        }
        // Report
        if (report)
            System.out.printf("%s %s %f : CRLB = %s, Devaitions = %s\n", solver.getClass().getSimpleName(), noiseModel, s, Arrays.toString(crlb), Arrays.toString(m.getStandardDeviation()));
    }
}