Examples with RandomDataGenerator org.apache.commons.math3.random.RandomDataGenerator used on opensource projects

Example 16 with RandomDataGenerator

use of org.apache.commons.math3.random.RandomDataGenerator in project GDSC-SMLM by aherbert.

the class PoissonCalculatorTest method canComputeLogLikelihoodRatio.

private void canComputeLogLikelihoodRatio(BaseNonLinearFunction nlf) {
    System.out.println(nlf.name);
    int n = maxx * maxx;
    double[] a = new double[] { 1 };
    // Simulate Poisson process
    nlf.initialise(a);
    RandomDataGenerator rdg = new RandomDataGenerator(new Well19937c(30051977));
    double[] x = new double[n];
    double[] u = new double[n];
    for (int i = 0; i < n; i++) {
        u[i] = nlf.eval(i);
        if (u[i] > 0)
            x[i] = rdg.nextPoisson(u[i]);
    }
    double ll = PoissonCalculator.logLikelihood(u, x);
    double mll = PoissonCalculator.maximumLogLikelihood(x);
    double llr = -2 * (ll - mll);
    double llr2 = PoissonCalculator.logLikelihoodRatio(u, x);
    System.out.printf("llr=%f, llr2=%f\n", llr, llr2);
    Assert.assertEquals("Log-likelihood ratio", llr, llr2, llr * 1e-10);
    double[] op = new double[x.length];
    for (int i = 0; i < n; i++) op[i] = PoissonCalculator.maximumLikelihood(x[i]);
    double max = Double.NEGATIVE_INFINITY;
    double maxa = 0;
    //TestSettings.setLogLevel(gdsc.smlm.TestSettings.LogLevel.DEBUG);
    int df = n - 1;
    ChiSquaredDistributionTable table = ChiSquaredDistributionTable.createUpperTailed(0.05, df);
    ChiSquaredDistributionTable table2 = ChiSquaredDistributionTable.createUpperTailed(0.001, df);
    System.out.printf("Chi2 = %f (q=%.3f), %f (q=%.3f)  %f %b  %f\n", table.getCrititalValue(df), table.getSignificanceValue(), table2.getCrititalValue(df), table2.getSignificanceValue(), ChiSquaredDistributionTable.computeQValue(24, 2), ChiSquaredDistributionTable.createUpperTailed(0.05, 2).reject(24, 2), ChiSquaredDistributionTable.getChiSquared(1e-6, 2));
    for (int i = 5; i <= 15; i++) {
        a[0] = (double) i / 10;
        nlf.initialise(a);
        for (int j = 0; j < n; j++) u[j] = nlf.eval(j);
        ll = PoissonCalculator.logLikelihood(u, x);
        llr = PoissonCalculator.logLikelihoodRatio(u, x);
        BigDecimal product = new BigDecimal(1);
        double ll2 = 0;
        for (int j = 0; j < n; j++) {
            double p1 = PoissonCalculator.likelihood(u[j], x[j]);
            ll2 += Math.log(p1);
            double ratio = p1 / op[j];
            product = product.multiply(new BigDecimal(ratio));
        }
        llr2 = -2 * Math.log(product.doubleValue());
        double p = ChiSquaredDistributionTable.computePValue(llr, df);
        double q = ChiSquaredDistributionTable.computeQValue(llr, df);
        TestSettings.info("a=%f, ll=%f, ll2=%f, llr=%f, llr2=%f, product=%s, p=%f, q=%f (reject=%b @ %.3f, reject=%b @ %.3f)\n", a[0], ll, ll2, llr, llr2, product.round(new MathContext(4)).toString(), p, q, table.reject(llr, df), table.getSignificanceValue(), table2.reject(llr, df), table2.getSignificanceValue());
        if (max < ll) {
            max = ll;
            maxa = a[0];
        }
        // too small to store in a double.
        if (product.doubleValue() > 0) {
            Assert.assertEquals("Log-likelihood", ll, ll2, Math.abs(ll2) * 1e-10);
            Assert.assertEquals("Log-likelihood ratio", llr, llr2, Math.abs(llr) * 1e-10);
        }
    }
    Assert.assertEquals("max", 1, maxa, 0);
}

Example 17 with RandomDataGenerator

use of org.apache.commons.math3.random.RandomDataGenerator in project GDSC-SMLM by aherbert.

the class PoissonGradientProcedureTest method gradientProcedureIsFasterUnrolledThanGradientProcedure.

private void gradientProcedureIsFasterUnrolledThanGradientProcedure(final int nparams, final boolean precomputed) {
    org.junit.Assume.assumeTrue(speedTests || TestSettings.RUN_SPEED_TESTS);
    final int iter = 100;
    rdg = new RandomDataGenerator(new Well19937c(30051977));
    final ArrayList<double[]> paramsList = new ArrayList<double[]>(iter);
    createFakeParams(nparams, iter, paramsList);
    // Remove the timing of the function call by creating a dummy function
    FakeGradientFunction f = new FakeGradientFunction(blockWidth, nparams);
    final Gradient1Function func = (precomputed) ? PrecomputedGradient1Function.wrapGradient1Function(f, Utils.newArray(f.size(), 0.1, 1.3)) : f;
    for (int i = 0; i < paramsList.size(); i++) {
        PoissonGradientProcedure p1 = new PoissonGradientProcedure(func);
        p1.computeFisherInformation(paramsList.get(i));
        p1.computeFisherInformation(paramsList.get(i));
        PoissonGradientProcedure p2 = PoissonGradientProcedureFactory.create(func);
        p2.computeFisherInformation(paramsList.get(i));
        p2.computeFisherInformation(paramsList.get(i));
        // Check they are the same
        Assert.assertArrayEquals("M " + i, p1.getLinear(), p2.getLinear(), 0);
    }
    // Realistic loops for an optimisation
    final int loops = 15;
    // Run till stable timing
    Timer t1 = new Timer() {

        @Override
        void run() {
            for (int i = 0, k = 0; i < paramsList.size(); i++) {
                PoissonGradientProcedure p1 = new PoissonGradientProcedure(func);
                for (int j = loops; j-- > 0; ) p1.computeFisherInformation(paramsList.get(k++ % iter));
            }
        }
    };
    long time1 = t1.getTime();
    Timer t2 = new Timer(t1.loops) {

        @Override
        void run() {
            for (int i = 0, k = 0; i < paramsList.size(); i++) {
                PoissonGradientProcedure p2 = PoissonGradientProcedureFactory.create(func);
                for (int j = loops; j-- > 0; ) p2.computeFisherInformation(paramsList.get(k++ % iter));
            }
        }
    };
    long time2 = t2.getTime();
    log("Precomputed=%b : Standard %d : Unrolled %d = %d : %fx\n", precomputed, time1, nparams, time2, (1.0 * time1) / time2);
    Assert.assertTrue(time2 < time1);
}

Example 18 with RandomDataGenerator

use of org.apache.commons.math3.random.RandomDataGenerator in project GDSC-SMLM by aherbert.

the class ChiSquaredDistributionTableTest method canPerformChiSquaredTest.

@Test
public void canPerformChiSquaredTest() {
    ChiSquareTest test = new ChiSquareTest();
    for (int n : new int[] { 10, 50, 100 }) {
        double[] x = Utils.newArray(n, 0.5, 1.0);
        long[] l = new long[x.length];
        RandomDataGenerator rdg = new RandomDataGenerator(new Well19937c(30051977));
        for (int i = 0; i < x.length; i++) l[i] = rdg.nextPoisson(x[i]);
        double chi2 = test.chiSquare(x, l);
        double ep = test.chiSquareTest(x, l);
        int df = x.length - 1;
        double o = ChiSquaredDistributionTable.computeQValue(chi2, df);
        Assert.assertEquals(ep, o, 1e-10);
        ChiSquaredDistributionTable upperTable = ChiSquaredDistributionTable.createUpperTailed(o, df);
        double upper = chi2 * 1.01;
        double lower = chi2 * 0.99;
        Assert.assertTrue("Upper did not reject higher", upperTable.reject(upper, df));
        Assert.assertFalse("Upper did not reject actual value", upperTable.reject(o, df));
        Assert.assertFalse("Upper did not accept lower", upperTable.reject(lower, df));
    }
}

Example 19 with RandomDataGenerator

use of org.apache.commons.math3.random.RandomDataGenerator in project GDSC-SMLM by aherbert.

the class FastMLEGradient2ProcedureTest method gradientProcedureComputesSameLogLikelihoodAsMLEGradientCalculator.

private void gradientProcedureComputesSameLogLikelihoodAsMLEGradientCalculator(int nparams) {
    int iter = 10;
    rdg = new RandomDataGenerator(new Well19937c(30051977));
    ArrayList<double[]> paramsList = new ArrayList<double[]>(iter);
    ArrayList<double[]> yList = new ArrayList<double[]>(iter);
    createFakeData(nparams, iter, paramsList, yList);
    FakeGradientFunction func = new FakeGradientFunction(blockWidth, nparams);
    MLEGradientCalculator calc = (MLEGradientCalculator) GradientCalculatorFactory.newCalculator(nparams, true);
    String name = String.format("[%d]", nparams);
    for (int i = 0; i < paramsList.size(); i++) {
        FastMLEGradient2Procedure p = FastMLEGradient2ProcedureFactory.createUnrolled(yList.get(i), func);
        double s = p.computeLogLikelihood(paramsList.get(i));
        double s2 = calc.logLikelihood(yList.get(i), paramsList.get(i), func);
        // Virtually the same ...
        Assert.assertEquals(name + " Result: Not same @ " + i, s, s2, Math.abs(s) * 1e-5);
    }
}

Example 20 with RandomDataGenerator

use of org.apache.commons.math3.random.RandomDataGenerator in project GDSC-SMLM by aherbert.

the class FastMLEGradient2ProcedureTest method gradientProcedureUnrolledComputesSameAsGradientProcedure.

private void gradientProcedureUnrolledComputesSameAsGradientProcedure(int nparams) {
    int iter = 10;
    rdg = new RandomDataGenerator(new Well19937c(30051977));
    ArrayList<double[]> paramsList = new ArrayList<double[]>(iter);
    ArrayList<double[]> yList = new ArrayList<double[]>(iter);
    createFakeData(nparams, iter, paramsList, yList);
    FakeGradientFunction func = new FakeGradientFunction(blockWidth, nparams);
    FastMLEGradient2Procedure p1, p2;
    String name = String.format("[%d]", nparams);
    for (int i = 0; i < paramsList.size(); i++) {
        p1 = new FastMLEGradient2Procedure(yList.get(i), func);
        p2 = FastMLEGradient2ProcedureFactory.createUnrolled(yList.get(i), func);
        double[] a = paramsList.get(i);
        double ll1 = p1.computeLogLikelihood(a);
        double ll2 = p2.computeLogLikelihood(a);
        Assert.assertEquals(name + " LL: Not same @ " + i, ll1, ll2, 0);
        p1 = new FastMLEGradient2Procedure(yList.get(i), func);
        p2 = FastMLEGradient2ProcedureFactory.createUnrolled(yList.get(i), func);
        p1.computeFirstDerivative(a);
        p2.computeFirstDerivative(a);
        Assert.assertArrayEquals(name + " first derivative value: Not same @ " + i, p1.u, p2.u, 0);
        Assert.assertArrayEquals(name + " first derivative: Not same @ " + i, p1.d1, p2.d1, 0);
        p1 = new FastMLEGradient2Procedure(yList.get(i), func);
        p2 = FastMLEGradient2ProcedureFactory.createUnrolled(yList.get(i), func);
        p1.computeSecondDerivative(a);
        p2.computeSecondDerivative(a);
        Assert.assertArrayEquals(name + " update value: Not same @ " + i, p1.u, p2.u, 0);
        Assert.assertArrayEquals(name + " update: Not same d1 @ " + i, p1.d1, p2.d1, 0);
        Assert.assertArrayEquals(name + " update: Not same d2 @ " + i, p1.d2, p2.d2, 0);
    }
}