use of uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.
the class ScmosLikelihoodWrapperTest method canComputePValue.
private static void canComputePValue(RandomSeed seed, BaseNonLinearFunction nlf) {
logger.log(TestLogUtils.getRecord(Level.INFO, nlf.name));
final int n = maxx * maxx;
final double[] a = new double[] { 1 };
// Simulate sCMOS camera
nlf.initialise(a);
final SCcmosLikelihoodWrapperTestData testData = (SCcmosLikelihoodWrapperTestData) dataCache.computeIfAbsent(seed, ScmosLikelihoodWrapperTest::createData);
final float[] var = testData.var;
final float[] g = testData.gain;
final float[] o = testData.offset;
final float[] sd = testData.sd;
final UniformRandomProvider r = RngUtils.create(seed.getSeed());
final SharedStateContinuousSampler gs = SamplerUtils.createGaussianSampler(r, 0, 1);
final double[] k = SimpleArrayUtils.newArray(n, 0, 1.0);
for (int i = 0; i < n; i++) {
double mean = nlf.eval(i);
if (mean > 0) {
mean = GdscSmlmTestUtils.createPoissonSampler(r, mean).sample();
}
k[i] = mean * g[i] + o[i] + gs.sample() * sd[i];
}
final ScmosLikelihoodWrapper f = new ScmosLikelihoodWrapper(nlf, a, k, n, var, g, o);
final double oll = f.computeObservedLikelihood();
double oll2 = 0;
final double[] op = new double[n];
for (int j = 0; j < n; j++) {
op[j] = ScmosLikelihoodWrapper.likelihood((k[j] - o[j]) / g[j], var[j], g[j], o[j], k[j]);
oll2 -= Math.log(op[j]);
}
logger.log(TestLogUtils.getRecord(Level.INFO, "oll=%f, oll2=%f", oll, oll2));
final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-10, 0);
TestAssertions.assertTest(oll2, oll, predicate, "Observed Log-likelihood");
final TDoubleArrayList list = new TDoubleArrayList();
final int imin = 5;
final int imax = 15;
for (int i = imin; i <= imax; i++) {
a[0] = (double) i / 10;
final double ll = f.likelihood(a);
list.add(ll);
final double llr = f.computeLogLikelihoodRatio(ll);
BigDecimal product = new BigDecimal(1);
double ll2 = 0;
for (int j = 0; j < n; j++) {
final double p1 = ScmosLikelihoodWrapper.likelihood(nlf.eval(j), var[j], g[j], o[j], k[j]);
ll2 -= Math.log(p1);
final double ratio = p1 / op[j];
product = product.multiply(new BigDecimal(ratio));
}
final double llr2 = -2 * Math.log(product.doubleValue());
final double q = f.computeQValue(ll);
logger.log(TestLogUtils.getRecord(Level.INFO, "a=%f, ll=%f, ll2=%f, llr=%f, llr2=%f, product=%s, p=%f", a[0], ll, ll2, llr, llr2, product.round(new MathContext(4)).toString(), q));
// too small to store in a double.
if (product.doubleValue() > 0) {
TestAssertions.assertTest(llr, llr2, predicate, "Log-likelihood");
}
}
// Find min using quadratic fit
final double[] data = list.toArray();
int index = SimpleArrayUtils.findMinIndex(data);
final double mina = (double) (imin + index) / 10;
double fita = mina;
try {
if (index == 0) {
index++;
}
if (index == data.length - 1) {
index--;
}
final int i1 = index - 1;
final int i2 = index;
final int i3 = index + 1;
fita = QuadraticUtils.findMinMax((double) (imin + i1) / 10, data[i1], (double) (imin + i2) / 10, data[i2], (double) (imin + i3) / 10, data[i3]);
} catch (final DataException ex) {
// Ignore
}
// Allow a tolerance as the random data may alter the p-value computation.
// Should allow it to be less than 2 increment either side of the answer.
logger.log(TestLogUtils.getRecord(Level.INFO, "min fit = %g => %g", mina, fita));
Assertions.assertEquals(1, fita, 0.199, "min");
}
use of uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.
the class PoissonGaussianConvolutionFunctionTest method probabilityMatchesLogProbability.
private static void probabilityMatchesLogProbability(final double gain, double mu, final double sd, boolean computePmf) {
// Note: The input s parameter is pre-gain.
final PoissonGaussianConvolutionFunction f = PoissonGaussianConvolutionFunction.createWithStandardDeviation(1.0 / gain, sd * gain);
f.setComputePmf(computePmf);
// Evaluate an initial range.
// Gaussian should have >99% within +/- s
// Poisson will have mean mu with a variance mu.
// At large mu it is approximately normal so use 3 sqrt(mu) for the range added to the mean
final int[] range = PoissonGaussianFunctionTest.getRange(gain, mu, sd);
final int min = range[0];
final int max = range[1];
// Note: The input mu parameter is pre-gain.
final double e = mu;
final Supplier<String> msg = () -> String.format("g=%f, mu=%f, s=%f, erf=%b", gain, mu, sd, computePmf);
final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-3, 0);
for (int x = min; x <= max; x++) {
final double p = f.likelihood(x, e);
if (p == 0) {
continue;
}
final double logP = f.logLikelihood(x, e);
TestAssertions.assertTest(Math.log(p), logP, predicate, msg);
}
}
use of uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.
the class FactorialTest method testFactorialDouble.
/**
* Test the factorial of a fractional number against Commons Math gamma(1+n).
*/
@SeededTest
void testFactorialDouble(RandomSeed seed) {
final UniformRandomProvider rng = RngUtils.create(seed.getSeed());
final DoubleDoubleBiPredicate tol = TestHelper.doublesAreClose(5e-15).or(TestHelper.doublesEqual());
for (int i = 0; i < 100; i++) {
final double n = rng.nextDouble() * 180;
final double expected = n < 1.5 ? 1 / (1 + Gamma.invGamma1pm1(n)) : Gamma.gamma(1 + n);
TestAssertions.assertTest(expected, Factorial.value(n), tol, () -> Double.toString(n));
}
}
use of uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.
the class FactorialTest method testFactorial.
/**
* Test the factorial of a fractional number against reference values.
*/
@ParameterizedTest
@CsvSource({ // @formatter:off
"0.25, 0.906402477055477077939", "0.75, 0.919062526848883233914", "1.125, 1.05946053733091417382", "1.75, 1.60835942198554565977", "2.75, 4.42298841046025056293", "3.25, 8.28508514183522016498", "4.75, 78.7844810613232131788", "6.25, 1155.3810139199896867", "12.75, 3255990905.16494153091", "17.125, 508919418354999.991547", "27.75, 1.32099626069721824877e+29", "57.75, 8.50381139447823243123e+77", "97.75, 2.99327649630298942593e+153", "137.75, 2.01698432357024766181e+236", "154.25, 1.08954758738660925068e+272", "164.25, 1.17747769003558155385e+294", "170.25, 2.62296687867878940345e+307", "170.624, 1.79421175992481035917e+308", // Limit
"170.6243769563027, 1.7976931348622301e+308", // Infinite
"170.62437695630274, 1.79769313486249261288e+308", "170.625, 1.80346206550076047216e+308", "171.5, 1.62639753771045308624e+310" // @formatter:on
})
void testFactorial(double n, double expected) {
final DoubleDoubleBiPredicate tol = TestHelper.doublesAreClose(1e-15).or(TestHelper.doublesEqual());
TestAssertions.assertTest(expected, Factorial.value(n), tol, () -> Double.toString(n));
}
use of uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.
the class MultivariateGaussianMixtureExpectationMaximizationTest method canCreateMixedMultivariateGaussianDistribution.
@SeededTest
void canCreateMixedMultivariateGaussianDistribution(RandomSeed seed) {
// Will be normalised
final double[] weights = { 1, 1 };
final double[][] means = new double[2][];
final double[][][] covariances = new double[2][][];
final double[][] data1 = { { 1, 2 }, { 2.5, 1.5 }, { 3.5, 1.0 } };
final double[][] data2 = { { 4, 2 }, { 3.5, -1.5 }, { -3.5, 1.0 } };
means[0] = getColumnMeans(data1);
covariances[0] = getCovariance(data1);
means[1] = getColumnMeans(data2);
covariances[1] = getCovariance(data2);
// Create components. This does not have to be zero based.
final LocalList<double[]> list = new LocalList<>();
list.addAll(Arrays.asList(data1));
list.addAll(Arrays.asList(data2));
final double[][] data = list.toArray(new double[0][]);
final int[] components = { -1, -1, -1, 3, 3, 3 };
// Randomise the data
for (int n = 0; n < 3; n++) {
final long start = n + seed.getSeedAsLong();
// This relies on the shuffle being the same
RandomUtils.shuffle(data, RngUtils.create(start));
RandomUtils.shuffle(components, RngUtils.create(start));
final MixtureMultivariateGaussianDistribution dist = MultivariateGaussianMixtureExpectationMaximization.createMixed(data, components);
Assertions.assertArrayEquals(new double[] { 0.5, 0.5 }, dist.getWeights());
final MultivariateGaussianDistribution[] distributions = dist.getDistributions();
Assertions.assertEquals(weights.length, distributions.length);
final DoubleDoubleBiPredicate test = TestHelper.doublesAreClose(1e-8);
for (int i = 0; i < means.length; i++) {
TestAssertions.assertArrayTest(means[i], distributions[i].getMeans(), test);
TestAssertions.assertArrayTest(covariances[i], distributions[i].getCovariances(), test);
}
}
}
Aggregations