Examples with CustomGammaDistribution org.apache.commons.math3.distribution.CustomGammaDistribution used on opensource projects

Example 1 with CustomGammaDistribution

use of org.apache.commons.math3.distribution.CustomGammaDistribution in project GDSC-SMLM by aherbert.

the class BaseFunctionSolverTest method drawGaussian.

/**
	 * Draw a Gaussian with Poisson shot noise and Gaussian read noise.
	 *
	 * @param params
	 *            The Gaussian parameters
	 * @param noise
	 *            The read noise
	 * @param noiseModel
	 *            the noise model
	 * @return The data
	 */
double[] drawGaussian(double[] params, double[] noise, NoiseModel noiseModel) {
    double[] data = new double[size * size];
    int n = params.length / 6;
    Gaussian2DFunction f = GaussianFunctionFactory.create2D(n, size, size, flags, null);
    f.initialise(params);
    // Poisson noise
    for (int i = 0; i < data.length; i++) {
        CustomPoissonDistribution dist = new CustomPoissonDistribution(randomGenerator, 1);
        double e = f.eval(i);
        if (e > 0) {
            dist.setMeanUnsafe(e);
            data[i] = dist.sample();
        }
    }
    // Simulate EM-gain
    if (noiseModel == NoiseModel.EMCCD) {
        // Use a gamma distribution
        // Since the call random.nextGamma(...) creates a Gamma distribution 
        // which pre-calculates factors only using the scale parameter we 
        // create a custom gamma distribution where the shape can be set as a property.
        CustomGammaDistribution dist = new CustomGammaDistribution(randomGenerator, 1, emGain);
        for (int i = 0; i < data.length; i++) {
            if (data[i] > 0) {
                dist.setShapeUnsafe(data[i]);
                // The sample will amplify the signal so we remap to the original scale
                data[i] = dist.sample() / emGain;
            }
        }
    }
    // Read-noise
    if (noise != null) {
        for (int i = 0; i < data.length; i++) {
            data[i] += randomGenerator.nextGaussian() * noise[i];
        }
    }
    //gdsc.core.ij.Utils.display("Spot", data, size, size);
    return data;
}

Example 2 with CustomGammaDistribution

use of org.apache.commons.math3.distribution.CustomGammaDistribution in project GDSC-SMLM by aherbert.

the class EMGainAnalysis method simulateFromPoissonGammaGaussian.

/**
	 * Randomly generate a histogram from poisson-gamma-gaussian samples
	 * 
	 * @return The histogram
	 */
private int[] simulateFromPoissonGammaGaussian() {
    // Randomly sample
    RandomGenerator random = new Well44497b(System.currentTimeMillis() + System.identityHashCode(this));
    PoissonDistribution poisson = new PoissonDistribution(random, _photons, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);
    CustomGammaDistribution gamma = new CustomGammaDistribution(random, _photons, _gain, GammaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
    final int steps = simulationSize;
    int[] sample = new int[steps];
    for (int n = 0; n < steps; n++) {
        if (n % 64 == 0)
            IJ.showProgress(n, steps);
        // Poisson
        double d = poisson.sample();
        // Gamma
        if (d > 0) {
            gamma.setShapeUnsafe(d);
            d = gamma.sample();
        }
        // Gaussian
        d += _noise * random.nextGaussian();
        // Convert the sample to a count 
        sample[n] = (int) Math.round(d + _bias);
    }
    int max = Maths.max(sample);
    int[] h = new int[max + 1];
    for (int s : sample) h[s]++;
    return h;
}