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));
}
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);
}
}
Aggregations