Examples with FisherInformationMatrix uk.ac.sussex.gdsc.smlm.fitting.FisherInformationMatrix used on opensource projects

Example 11 with FisherInformationMatrix

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

the class PoissonGradientProcedureTest method varianceMatchesFormula.

@SeededTest
void varianceMatchesFormula() {
    // Assumptions.assumeTrue(false);
    final double[] n_ = new double[] { 20, 50, 100, 500 };
    final double[] b2_ = new double[] { 0, 1, 2, 4 };
    final double[] s_ = new double[] { 1, 1.2, 1.5 };
    final double[] x_ = new double[] { 4.8, 5, 5.5 };
    final double a = 100;
    final int size = 10;
    final Gaussian2DFunction f = GaussianFunctionFactory.create2D(1, size, size, GaussianFunctionFactory.FIT_ERF_CIRCLE, null);
    final PoissonGradientProcedure p = PoissonGradientProcedureUtils.create(f);
    final int ix = f.findGradientIndex(Gaussian2DFunction.X_POSITION);
    final int iy = f.findGradientIndex(Gaussian2DFunction.Y_POSITION);
    final double[] params = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    for (final double n : n_) {
        params[Gaussian2DFunction.SIGNAL] = n;
        for (final double b2 : b2_) {
            params[Gaussian2DFunction.BACKGROUND] = b2;
            for (final double s : s_) {
                final double ss = s * a;
                params[Gaussian2DFunction.X_SD] = s;
                for (final double x : x_) {
                    params[Gaussian2DFunction.X_POSITION] = x;
                    for (final double y : x_) {
                        params[Gaussian2DFunction.Y_POSITION] = y;
                        p.computeFisherInformation(params);
                        final FisherInformationMatrix m1 = new FisherInformationMatrix(p.getLinear(), p.numberOfGradients);
                        final double[] crlb = m1.crlb();
                        if (crlb == null) {
                            Assertions.fail("No variance");
                        }
                        @SuppressWarnings("null") final double o1 = Math.sqrt(crlb[ix]) * a;
                        final double o2 = Math.sqrt(crlb[iy]) * a;
                        final double e = Gaussian2DPeakResultHelper.getMLPrecisionX(a, ss, n, b2, false);
                        // logger.fine(FunctionUtils.getSupplier("e = %f : o = %f %f", e, o1, o2);
                        Assertions.assertEquals(e, o1, e * 5e-2);
                        Assertions.assertEquals(e, o2, e * 5e-2);
                    }
                }
            }
        }
    }
}

Example 12 with FisherInformationMatrix

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

the class BaseFunctionSolverTest method canFitSingleGaussian.

void canFitSingleGaussian(RandomSeed seed, FunctionSolver solver, boolean applyBounds, NoiseModel noiseModel) {
    // Allow reporting the fit deviations
    final boolean report = false;
    double[] crlb = null;
    SimpleArrayMoment moment = null;
    final double[] noise = getNoise(seed, noiseModel);
    if (solver.isWeighted()) {
        solver.setWeights(getWeights(seed, noiseModel));
    }
    final UniformRandomProvider rg = RngUtils.create(seed.getSeed());
    for (final double s : signal) {
        final double[] expected = createParams(1, s, 0, 0, 1);
        final double[] lower = createParams(0, s * 0.5, -0.3, -0.3, 0.8);
        final double[] upper = createParams(3, s * 2, 0.3, 0.3, 1.2);
        if (applyBounds) {
            solver.setBounds(lower, upper);
        }
        if (report) {
            // Compute the CRLB for a Poisson process
            final PoissonGradientProcedure gp = PoissonGradientProcedureUtils.create((Gradient1Function) ((BaseFunctionSolver) solver).getGradientFunction());
            gp.computeFisherInformation(expected);
            final FisherInformationMatrix f = new FisherInformationMatrix(gp.getLinear(), gp.numberOfGradients);
            crlb = f.crlbSqrt();
            // Compute the deviations.
            // Note this is not the same as the CRLB as the fit is repeated
            // against the same data.
            // It should be repeated against different data generated with constant
            // parameters and variable noise.
            moment = new SimpleArrayMoment();
        }
        final double[] data = drawGaussian(expected, noise, noiseModel, rg);
        for (final double db : base) {
            for (final double dx : shift) {
                for (final double dy : shift) {
                    for (final double dsx : factor) {
                        final double[] p = createParams(db, s, dx, dy, dsx);
                        final double[] fp = fitGaussian(solver, data, p, expected);
                        for (int i = 0; i < expected.length; i++) {
                            if (fp[i] < lower[i]) {
                                Assertions.fail(FunctionUtils.getSupplier("Fit Failed: [%d] %.2f < %.2f: %s != %s", i, fp[i], lower[i], Arrays.toString(fp), Arrays.toString(expected)));
                            }
                            if (fp[i] > upper[i]) {
                                Assertions.fail(FunctionUtils.getSupplier("Fit Failed: [%d] %.2f > %.2f: %s != %s", i, fp[i], upper[i], Arrays.toString(fp), Arrays.toString(expected)));
                            }
                            if (report) {
                                fp[i] = expected[i] - fp[i];
                            }
                        }
                        // Store the deviations
                        if (report) {
                            moment.add(fp);
                        }
                    }
                }
            }
        }
        // Report
        if (report) {
            logger.info(FunctionUtils.getSupplier("%s %s %f : CRLB = %s, Deviations = %s", solver.getClass().getSimpleName(), noiseModel, s, Arrays.toString(crlb), Arrays.toString(moment.getStandardDeviation())));
        }
    }
}

Example 13 with FisherInformationMatrix

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

the class MleLvmSteppingFunctionSolver method computeLastFisherInformationMatrix.

@Override
protected FisherInformationMatrix computeLastFisherInformationMatrix(double[] fx) {
    // The Hessian matrix refers to the log-likelihood ratio.
    // Compute and invert a matrix related to the Poisson log-likelihood.
    // This assumes this does achieve the maximum likelihood estimate for a
    // Poisson process.
    Gradient1Function localF1 = (Gradient1Function) function;
    // Capture the y-values if necessary
    if (fx != null && fx.length == localF1.size()) {
        localF1 = new Gradient2FunctionValueStore(localF1, fx);
    }
    // Add the weights if necessary
    if (weights != null) {
        localF1 = OffsetGradient1Function.wrapGradient1Function(localF1, weights);
    }
    final PoissonGradientProcedure p = PoissonGradientProcedureUtils.create(localF1);
    p.computeFisherInformation(lastA);
    if (p.isNaNGradients()) {
        throw new FunctionSolverException(FitStatus.INVALID_GRADIENTS);
    }
    // Re-use space
    p.getLinear(walpha);
    return new FisherInformationMatrix(walpha, beta.length);
}

Example 14 with FisherInformationMatrix

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

the class MleLvmSteppingFunctionSolver method computeFunctionFisherInformationMatrix.

@Override
protected FisherInformationMatrix computeFunctionFisherInformationMatrix(double[] y, double[] a) {
    // Compute and invert a matrix related to the Poisson log-likelihood.
    // This assumes this does achieve the maximum likelihood estimate for a
    // Poisson process.
    // We must wrap the gradient function if weights are present.
    Gradient1Function localF1 = (Gradient1Function) function;
    if (weights != null) {
        localF1 = OffsetGradient1Function.wrapGradient1Function(localF1, weights);
    }
    final PoissonGradientProcedure p = PoissonGradientProcedureUtils.create(localF1);
    p.computeFisherInformation(a);
    if (p.isNaNGradients()) {
        throw new FunctionSolverException(FitStatus.INVALID_GRADIENTS);
    }
    return new FisherInformationMatrix(p.getLinear(), function.getNumberOfGradients());
}

Example 15 with FisherInformationMatrix

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

the class NonLinearFit method computeFitDeviations.

/**
 * Compute the parameter deviations using the covariance matrix of the solution.
 *
 * @param n the n
 * @param parametersVariance the parameter deviations
 * @return true, if successful
 */
private boolean computeFitDeviations(int n, double[] parametersVariance) {
    if (isMle()) {
        // The Hessian matrix refers to the log-likelihood ratio.
        // Compute and invert a matrix related to the Poisson log-likelihood.
        // This assumes this does achieve the maximum likelihood estimate for a
        // Poisson process.
        final double[][] fim = calculator.fisherInformationMatrix(n, null, func);
        if (calculator.isNaNGradients()) {
            throw new FunctionSolverException(FitStatus.INVALID_GRADIENTS);
        }
        // Use a dedicated solver optimised for inverting the matrix diagonal.
        final FisherInformationMatrix m = new FisherInformationMatrix(fim);
        setDeviations(parametersVariance, m);
        return true;
    }
    final double[] var = calculator.variance(n, null, func);
    if (var != null) {
        setDeviations(parametersVariance, var);
        return true;
    }
    return false;
}