Search in sources :

Example 1 with FisherInformation

use of uk.ac.sussex.gdsc.smlm.function.FisherInformation in project GDSC-SMLM by aherbert.

the class UnivariateLikelihoodFisherInformationCalculatorTest method canComputePerPixelPoissonGaussianApproximationFisherInformation.

private static void canComputePerPixelPoissonGaussianApproximationFisherInformation(UniformRandomProvider rng) {
    // Create function
    final Gaussian2DFunction func = GaussianFunctionFactory.create2D(1, 10, 10, GaussianFunctionFactory.FIT_ERF_CIRCLE, null);
    final double[] params = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    params[Gaussian2DFunction.BACKGROUND] = nextUniform(rng, 0.1, 0.3);
    params[Gaussian2DFunction.SIGNAL] = nextUniform(rng, 100, 300);
    params[Gaussian2DFunction.X_POSITION] = nextUniform(rng, 4, 6);
    params[Gaussian2DFunction.Y_POSITION] = nextUniform(rng, 4, 6);
    params[Gaussian2DFunction.X_SD] = nextUniform(rng, 1, 1.3);
    Gradient1Function f1 = func;
    FisherInformation[] fi;
    // Get a per-pixel variance
    final double[] var = new double[func.size()];
    fi = new FisherInformation[var.length];
    for (int i = var.length; i-- > 0; ) {
        var[i] = 0.9 + 0.2 * rng.nextDouble();
        fi[i] = new PoissonGaussianApproximationFisherInformation(Math.sqrt(var[i]));
    }
    f1 = (Gradient1Function) OffsetFunctionFactory.wrapFunction(func, var);
    // This introduces a dependency on a different package, and relies on that
    // computing the correct answer. However that code predates this and so the
    // test ensures that the FisherInformationCalculator functions correctly.
    final PoissonGradientProcedure p1 = PoissonGradientProcedureUtils.create(f1);
    p1.computeFisherInformation(params);
    final double[] e = p1.getLinear();
    final FisherInformationCalculator calc = new UnivariateLikelihoodFisherInformationCalculator(func, fi);
    final FisherInformationMatrix I = calc.compute(params);
    final double[] o = I.getMatrix().data;
    TestAssertions.assertArrayTest(e, o, TestHelper.doublesAreClose(1e-6, 0));
}
Also used : Gradient1Function(uk.ac.sussex.gdsc.smlm.function.Gradient1Function) Gaussian2DFunction(uk.ac.sussex.gdsc.smlm.function.gaussian.Gaussian2DFunction) PoissonGradientProcedure(uk.ac.sussex.gdsc.smlm.fitting.nonlinear.gradient.PoissonGradientProcedure) PoissonFisherInformation(uk.ac.sussex.gdsc.smlm.function.PoissonFisherInformation) PoissonGaussianApproximationFisherInformation(uk.ac.sussex.gdsc.smlm.function.PoissonGaussianApproximationFisherInformation) HalfPoissonFisherInformation(uk.ac.sussex.gdsc.smlm.function.HalfPoissonFisherInformation) FisherInformation(uk.ac.sussex.gdsc.smlm.function.FisherInformation) PoissonGaussianApproximationFisherInformation(uk.ac.sussex.gdsc.smlm.function.PoissonGaussianApproximationFisherInformation)

Example 2 with FisherInformation

use of uk.ac.sussex.gdsc.smlm.function.FisherInformation in project GDSC-SMLM by aherbert.

the class UnivariateLikelihoodFisherInformationCalculatorTest method computePoissonFisherInformation.

private static void computePoissonFisherInformation(UniformRandomProvider rng, Model model) {
    // Create function
    final Gaussian2DFunction func = GaussianFunctionFactory.create2D(1, 10, 10, GaussianFunctionFactory.FIT_ERF_CIRCLE, null);
    final double[] params = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    params[Gaussian2DFunction.BACKGROUND] = nextUniform(rng, 0.1, 0.3);
    params[Gaussian2DFunction.SIGNAL] = nextUniform(rng, 100, 300);
    params[Gaussian2DFunction.X_POSITION] = nextUniform(rng, 4, 6);
    params[Gaussian2DFunction.Y_POSITION] = nextUniform(rng, 4, 6);
    params[Gaussian2DFunction.X_SD] = nextUniform(rng, 1, 1.3);
    Gradient1Function f1 = func;
    FisherInformation fi;
    switch(model) {
        // Get a variance
        case POISSON_GAUSSIAN:
            final double var = 0.9 + 0.2 * rng.nextDouble();
            fi = new PoissonGaussianApproximationFisherInformation(Math.sqrt(var));
            f1 = (Gradient1Function) OffsetFunctionFactory.wrapFunction(func, SimpleArrayUtils.newDoubleArray(func.size(), var));
            break;
        case POISSON:
            fi = new PoissonFisherInformation();
            break;
        case HALF_POISSON:
            fi = new HalfPoissonFisherInformation();
            break;
        default:
            throw new IllegalStateException();
    }
    // This introduces a dependency on a different package, and relies on that
    // computing the correct answer. However that code predates this and so the
    // test ensures that the FisherInformationCalculator functions correctly.
    final PoissonGradientProcedure p1 = PoissonGradientProcedureUtils.create(f1);
    p1.computeFisherInformation(params);
    final double[] e = p1.getLinear();
    final FisherInformationCalculator calc = new UnivariateLikelihoodFisherInformationCalculator(func, fi);
    final FisherInformationMatrix I = calc.compute(params);
    final double[] o = I.getMatrix().data;
    final boolean emCcd = model == Model.HALF_POISSON;
    if (emCcd) {
        // Assumes half the poisson fisher information
        SimpleArrayUtils.multiply(e, 0.5);
    }
    Assertions.assertArrayEquals(e, o, 1e-6);
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(5e-2, 0);
    if (model == Model.POISSON || model == Model.HALF_POISSON) {
        // Get the Mortensen approximation for fitting Poisson data with a Gaussian.
        // Set a to 100 for the square pixel adjustment.
        final double a = 100;
        final double s = params[Gaussian2DFunction.X_SD] * a;
        final double N = params[Gaussian2DFunction.SIGNAL];
        final double b2 = params[Gaussian2DFunction.BACKGROUND];
        double var = Gaussian2DPeakResultHelper.getMLVarianceX(a, s, N, b2, emCcd);
        // Convert expected variance to pixels
        var /= (a * a);
        // Get the limits by inverting the Fisher information
        final double[] crlb = I.crlb();
        TestAssertions.assertTest(var, crlb[2], predicate);
        TestAssertions.assertTest(var, crlb[3], predicate);
    }
}
Also used : DoubleDoubleBiPredicate(uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate) Gradient1Function(uk.ac.sussex.gdsc.smlm.function.Gradient1Function) PoissonGradientProcedure(uk.ac.sussex.gdsc.smlm.fitting.nonlinear.gradient.PoissonGradientProcedure) PoissonGaussianApproximationFisherInformation(uk.ac.sussex.gdsc.smlm.function.PoissonGaussianApproximationFisherInformation) HalfPoissonFisherInformation(uk.ac.sussex.gdsc.smlm.function.HalfPoissonFisherInformation) PoissonFisherInformation(uk.ac.sussex.gdsc.smlm.function.PoissonFisherInformation) HalfPoissonFisherInformation(uk.ac.sussex.gdsc.smlm.function.HalfPoissonFisherInformation) Gaussian2DFunction(uk.ac.sussex.gdsc.smlm.function.gaussian.Gaussian2DFunction) PoissonFisherInformation(uk.ac.sussex.gdsc.smlm.function.PoissonFisherInformation) PoissonGaussianApproximationFisherInformation(uk.ac.sussex.gdsc.smlm.function.PoissonGaussianApproximationFisherInformation) HalfPoissonFisherInformation(uk.ac.sussex.gdsc.smlm.function.HalfPoissonFisherInformation) FisherInformation(uk.ac.sussex.gdsc.smlm.function.FisherInformation)

Aggregations

PoissonGradientProcedure (uk.ac.sussex.gdsc.smlm.fitting.nonlinear.gradient.PoissonGradientProcedure)2 FisherInformation (uk.ac.sussex.gdsc.smlm.function.FisherInformation)2 Gradient1Function (uk.ac.sussex.gdsc.smlm.function.Gradient1Function)2 HalfPoissonFisherInformation (uk.ac.sussex.gdsc.smlm.function.HalfPoissonFisherInformation)2 PoissonFisherInformation (uk.ac.sussex.gdsc.smlm.function.PoissonFisherInformation)2 PoissonGaussianApproximationFisherInformation (uk.ac.sussex.gdsc.smlm.function.PoissonGaussianApproximationFisherInformation)2 Gaussian2DFunction (uk.ac.sussex.gdsc.smlm.function.gaussian.Gaussian2DFunction)2 DoubleDoubleBiPredicate (uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate)1