Java RandomDataGenerator.nextPoisson方法代碼示例

import org.apache.commons.math3.random.RandomDataGenerator; //導入方法依賴的package包/類
/**
 * Get the number of blinks using the specified random data generator using a Poisson or Geometric distribution.
 * @param useGeometricBlinkingDistribution
 * @param rand
 * @param mean
 * @return The number of blinks
 */
public static int getBlinks(boolean useGeometricBlinkingDistribution, RandomDataGenerator rand, double mean)
{
	if (mean > 0)
	{
		return (useGeometricBlinkingDistribution) ? nextGeometric(rand, mean) : (int) rand.nextPoisson(mean);
	}
	return 0;
}
 

import org.apache.commons.math3.random.RandomDataGenerator; //導入方法依賴的package包/類
@Test
public void mleCalculatorComputesLogLikelihoodRatio()
{
	EllipticalGaussian2DFunction func = new EllipticalGaussian2DFunction(1, blockWidth, blockWidth);
	int n = blockWidth * blockWidth;
	double[] a = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
	rdg = new RandomDataGenerator(new Well19937c(30051977));
	for (int run = 5; run-- > 0;)
	{
		a[Gaussian2DFunction.BACKGROUND] = random(Background);
		a[Gaussian2DFunction.SIGNAL] = random(Amplitude);
		a[Gaussian2DFunction.ANGLE] = random(Angle);
		a[Gaussian2DFunction.X_POSITION] = random(Xpos);
		a[Gaussian2DFunction.Y_POSITION] = random(Ypos);
		a[Gaussian2DFunction.X_SD] = random(Xwidth);
		a[Gaussian2DFunction.Y_SD] = random(Ywidth);

		// Simulate Poisson process
		func.initialise(a);
		double[] x = SimpleArrayUtils.newArray(n, 0, 1.0);
		double[] u = new double[x.length];
		for (int i = 0; i < n; i++)
		{
			u[i] = func.eval(i);
			if (u[i] > 0)
				x[i] = rdg.nextPoisson(u[i]);
		}

		int ng = func.getNumberOfGradients();
		double[][] alpha = new double[ng][ng];
		double[] beta = new double[ng];

		GradientCalculator calc = GradientCalculatorFactory.newCalculator(ng, true);

		double llr = PoissonCalculator.logLikelihoodRatio(u, x);
		double llr2 = calc.findLinearised(n, x, a, alpha, beta, func);
		//System.out.printf("llr=%f, llr2=%f\n", llr, llr2);
		Assert.assertEquals("Log-likelihood ratio", llr, llr2, llr * 1e-10);
	}
}
 

import org.apache.commons.math3.random.RandomDataGenerator; //導入方法依賴的package包/類
@Test
public void canPerformChiSquaredTest()
{
	ChiSquareTest test = new ChiSquareTest();
	for (int n : new int[] { 10, 50, 100 })
	{
		double[] x = SimpleArrayUtils.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));
	}
}
 

import org.apache.commons.math3.random.RandomDataGenerator; //導入方法依賴的package包/類
private static short[] combine(float[] data1, float[] data2, float[] data3)
{
	// Combine images and add a bias and read noise
	RandomDataGenerator rg = new RandomDataGenerator(rand);
	short[] data = new short[data1.length];
	for (int i = 0; i < data.length; i++)
	{
		final double mu = data1[i] + data2[i] + data3[i];
		data[i] = (short) (((mu != 0) ? rg.nextPoisson(mu) : 0) + rg.nextGaussian(bias, 5));
	}
	return data;
}
 

import org.apache.commons.math3.random.RandomDataGenerator; //導入方法依賴的package包/類
private void cannotSubtractConstantBackgroundAndComputeLogLikelihoodRatio(BaseNonLinearFunction nlf1,
		BaseNonLinearFunction nlf2, BaseNonLinearFunction nlf3)
{
	//System.out.println(nlf1.name);

	int n = maxx * maxx;

	double[] a = new double[] { 1 };

	// Simulate Poisson process of 3 combined functions
	nlf1.initialise(a);
	nlf2.initialise(a);
	nlf3.initialise(a);
	RandomDataGenerator rdg = new RandomDataGenerator(new Well19937c(30051977));
	double[] x = SimpleArrayUtils.newArray(n, 0, 1.0);
	double[] u = new double[x.length];
	double[] b1 = new double[x.length];
	double[] b2 = new double[x.length];
	double[] b3 = new double[x.length];
	for (int i = 0; i < n; i++)
	{
		b1[i] = nlf1.eval(i);
		b2[i] = nlf2.eval(i);
		b3[i] = nlf3.eval(i);
		u[i] = b1[i] + b2[i] + b3[i];
		if (u[i] > 0)
			x[i] = rdg.nextPoisson(u[i]);
	}

	// x is the target data
	// b1 is the already computed background
	// b2 is the first function to add to the model
	// b3 is the second function to add to the model

	// Compute the LLR of adding b3 to b2 given we already have b1 modelling data x
	double[] b12 = add(b1, b2);
	double ll1a = PoissonCalculator.logLikelihood(b12, x);
	double ll2a = PoissonCalculator.logLikelihood(add(b12, b3), x);
	double llra = -2 * (ll1a - ll2a);
	//System.out.printf("x|(a+b+c) ll1=%f, ll2=%f, llra=%f\n", ll1a, ll2a, llra);

	// Compute the LLR of adding b3 to b2 given we already have x minus b1
	x = subtract(x, b1);
	double ll1b = PoissonCalculator.logLikelihood(b2, x);
	double ll2b = PoissonCalculator.logLikelihood(add(b2, b3), x);
	double llrb = -2 * (ll1b - ll2b);
	//System.out.printf("x-a|(b+c) : ll1=%f, ll2=%f, llrb=%f\n", ll1b, ll2b, llrb);

	//System.out.printf("llr=%f (%g), llr2=%f (%g)\n", llra, PoissonCalculator.computePValue(llra, 1), llrb,
	//		PoissonCalculator.computePValue(llrb, 1));
	Assert.assertNotEquals("Log-likelihood ratio", llra, llrb, llra * 1e-10);

}
 

import org.apache.commons.math3.random.RandomDataGenerator; //導入方法依賴的package包/類
private void functionComputesTargetGradient(Gaussian2DFunction f1, int targetParameter, double threshold)
{
	int[] indices = f1.gradientIndices();
	int gradientIndex = findGradientIndex(f1, targetParameter);
	double[] dyda = new double[indices.length];
	double[] a;

	SCMOSLikelihoodWrapper ff1;

	int n = maxx * maxx;
	int count = 0, total = 0;

	RandomDataGenerator rdg = new RandomDataGenerator(new Well19937c(30051977));

	for (double background : testbackground)
		for (double signal1 : testsignal1)
			for (double cx1 : testcx1)
				for (double cy1 : testcy1)
					for (double cz1 : testcz1)
						for (double[] w1 : testw1)
							for (double angle1 : testangle1)
							{
								a = createParameters(background, signal1, cx1, cy1, cz1, w1[0], w1[1], angle1);

								// Create y as a function we would want to move towards
								double[] a2 = a.clone();
								a2[targetParameter] *= 1.3;
								f1.initialise(a2);
								double[] data = new double[n];
								for (int i = 0; i < n; i++)
								{
									// Simulate sCMOS camera
									double u = f1.eval(i);
									data[i] = rdg.nextPoisson(u) * g[i] + rdg.nextGaussian(o[i], sd[i]);
								}

								ff1 = new SCMOSLikelihoodWrapper(f1, a, data, n, var, g, o);

								// Numerically solve gradient. 
								// Calculate the step size h to be an exact numerical representation
								final double xx = a[targetParameter];

								// Get h to minimise roundoff error
								double h = Precision.representableDelta(xx, h_);

								ff1.likelihood(getVariables(indices, a), dyda);

								// Evaluate at (x+h) and (x-h)
								a[targetParameter] = xx + h;
								double value2 = ff1.likelihood(getVariables(indices, a));

								a[targetParameter] = xx - h;
								double value3 = ff1.likelihood(getVariables(indices, a));

								double gradient = (value2 - value3) / (2 * h);
								boolean ok = Math.signum(gradient) == Math.signum(dyda[gradientIndex]) ||
										Math.abs(gradient - dyda[gradientIndex]) < 0.1;
								//logf("[%s-%s]/2*%g : %g == %g\n", "" + value2, "" + value3, h, gradient,
								//		dyda[gradientIndex]);
								if (!ok)
									Assert.assertTrue(
											NAME[targetParameter] + ": " + gradient + " != " + dyda[gradientIndex],
											ok);
								ok = eq.almostEqualRelativeOrAbsolute(gradient, dyda[gradientIndex]);
								if (ok)
									count++;
								total++;

							}
	double p = (100.0 * count) / total;
	logf("%s : %s = %d / %d (%.2f)\n", f1.getClass().getSimpleName(), NAME[targetParameter], count, total, p);
	Assert.assertTrue(NAME[targetParameter] + " fraction too low: " + p, p > threshold);
}
 

import org.apache.commons.math3.random.RandomDataGenerator; //導入方法依賴的package包/類
private void canComputePValue(BaseNonLinearFunction nlf)
{
	System.out.println(nlf.name);

	int n = maxx * maxx;

	double[] a = new double[] { 1 };

	// Simulate sCMOS camera
	nlf.initialise(a);
	RandomDataGenerator rdg = new RandomDataGenerator(new Well19937c(30051977));
	double[] k = SimpleArrayUtils.newArray(n, 0, 1.0);
	for (int i = 0; i < n; i++)
	{
		double u = nlf.eval(i);
		if (u > 0)
			u = rdg.nextPoisson(u);
		k[i] = u * g[i] + rdg.nextGaussian(o[i], sd[i]);
	}

	SCMOSLikelihoodWrapper f = new SCMOSLikelihoodWrapper(nlf, a, k, n, var, g, o);

	double oll = f.computeObservedLikelihood();
	double oll2 = 0;
	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]);
	}
	System.out.printf("oll=%f, oll2=%f\n", oll, oll2);
	Assert.assertEquals("Observed Log-likelihood", oll2, oll, oll2 * 1e-10);

	double min = Double.POSITIVE_INFINITY;
	double mina = 0;
	for (int i = 5; i <= 15; i++)
	{
		a[0] = (double) i / 10;
		double ll = f.likelihood(a);
		double llr = f.computeLogLikelihoodRatio(ll);
		BigDecimal product = new BigDecimal(1);
		double ll2 = 0;
		for (int j = 0; j < n; j++)
		{
			double p1 = SCMOSLikelihoodWrapper.likelihood(nlf.eval(j), var[j], g[j], o[j], k[j]);
			ll2 -= Math.log(p1);
			double ratio = p1 / op[j];
			product = product.multiply(new BigDecimal(ratio));
		}
		double llr2 = -2 * Math.log(product.doubleValue());
		double q = f.computeQValue(ll);
		System.out.printf("a=%f, ll=%f, ll2=%f, llr=%f, llr2=%f, product=%s, p=%f\n", a[0], ll, ll2, llr, llr2,
				product.round(new MathContext(4)).toString(), q);
		if (min > ll)
		{
			min = ll;
			mina = a[0];
		}

		// Only value if the product could be computed. Low ratios cause it to becomes 
		// too small to store in a double.
		if (product.doubleValue() > 0)
			Assert.assertEquals("Log-likelihood", llr, llr2, llr * 1e-10);
	}

	Assert.assertEquals("min", 1, mina, 0);
}
 

import org.apache.commons.math3.random.RandomDataGenerator; //導入方法依賴的package包/類
/**
 * Get random data from a Poisson distribution
 * @param mean Poisson mean
 * @return Random value
 */
public static double poisson(double mean){
    RandomDataGenerator rdg = new RandomDataGenerator();
    return rdg.nextPoisson(mean);
}