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

Example 1 with ExponentialDistribution

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

the class CMOSAnalysis method simulate.

private void simulate() {
    // Create the offset, variance and gain for each pixel
    int n = size * size;
    float[] pixelOffset = new float[n];
    float[] pixelVariance = new float[n];
    float[] pixelGain = new float[n];
    IJ.showStatus("Creating random per-pixel readout");
    long start = System.currentTimeMillis();
    RandomGenerator rg = new Well19937c();
    PoissonDistribution pd = new PoissonDistribution(rg, offset, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);
    ExponentialDistribution ed = new ExponentialDistribution(rg, variance, ExponentialDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
    totalProgress = n;
    stepProgress = Utils.getProgressInterval(totalProgress);
    for (int i = 0; i < n; i++) {
        if (i % n == 0)
            IJ.showProgress(i, n);
        // Q. Should these be clipped to a sensible range?
        pixelOffset[i] = (float) pd.sample();
        pixelVariance[i] = (float) ed.sample();
        pixelGain[i] = (float) (gain + rg.nextGaussian() * gainSD);
    }
    IJ.showProgress(1);
    // Avoid all the file saves from updating the progress bar and status line
    Utils.setShowStatus(false);
    Utils.setShowProgress(false);
    JLabel statusLine = Utils.getStatusLine();
    progressBar = Utils.getProgressBar();
    // Save to the directory as a stack
    ImageStack simulationStack = new ImageStack(size, size);
    simulationStack.addSlice("Offset", pixelOffset);
    simulationStack.addSlice("Variance", pixelVariance);
    simulationStack.addSlice("Gain", pixelGain);
    simulationImp = new ImagePlus("PerPixel", simulationStack);
    // Only the info property is saved to the TIFF file
    simulationImp.setProperty("Info", String.format("Offset=%s; Variance=%s; Gain=%s +/- %s", Utils.rounded(offset), Utils.rounded(variance), Utils.rounded(gain), Utils.rounded(gainSD)));
    IJ.save(simulationImp, new File(directory, "perPixelSimulation.tif").getPath());
    // Create thread pool and workers
    ExecutorService executor = Executors.newFixedThreadPool(getThreads());
    TurboList<Future<?>> futures = new TurboList<Future<?>>(nThreads);
    // Simulate the zero exposure input.
    // Simulate 20 - 200 photon images.
    int[] photons = new int[] { 0, 20, 50, 100, 200 };
    totalProgress = photons.length * frames;
    stepProgress = Utils.getProgressInterval(totalProgress);
    progress = 0;
    progressBar.show(0);
    // For saving stacks
    int blockSize = 10;
    int nPerThread = (int) Math.ceil((double) frames / nThreads);
    // Convert to fit the block size
    nPerThread = (int) Math.ceil((double) nPerThread / blockSize) * blockSize;
    long seed = start;
    for (int p : photons) {
        statusLine.setText("Simulating " + Utils.pleural(p, "photon"));
        // Create the directory
        File out = new File(directory, String.format("photon%03d", p));
        if (!out.exists())
            out.mkdir();
        for (int from = 0; from < frames; ) {
            int to = Math.min(from + nPerThread, frames);
            futures.add(executor.submit(new SimulationWorker(seed++, out.getPath(), simulationStack, from, to, blockSize, p)));
            from = to;
        }
        // Wait for all to finish
        for (int t = futures.size(); t-- > 0; ) {
            try {
                // The future .get() method will block until completed
                futures.get(t).get();
            } catch (Exception e) {
                // This should not happen. 
                e.printStackTrace();
            }
        }
        futures.clear();
    }
    Utils.setShowStatus(true);
    Utils.setShowProgress(true);
    IJ.showProgress(1);
    executor.shutdown();
    Utils.log("Simulation time = " + Utils.timeToString(System.currentTimeMillis() - start));
}

Example 2 with ExponentialDistribution

use of org.apache.commons.math3.distribution.ExponentialDistribution in project deeplearning4j by deeplearning4j.

the class TestReconstructionDistributions method testExponentialLogProb.

@Test
public void testExponentialLogProb() {
    Nd4j.getRandom().setSeed(12345);
    int inputSize = 4;
    int[] mbs = new int[] { 1, 2, 5 };
    Random r = new Random(12345);
    for (boolean average : new boolean[] { true, false }) {
        for (int minibatch : mbs) {
            INDArray x = Nd4j.zeros(minibatch, inputSize);
            for (int i = 0; i < minibatch; i++) {
                for (int j = 0; j < inputSize; j++) {
                    x.putScalar(i, j, r.nextInt(2));
                }
            }
            //i.e., pre-afn gamma
            INDArray distributionParams = Nd4j.rand(minibatch, inputSize).muli(2).subi(1);
            INDArray gammas = Transforms.tanh(distributionParams, true);
            ReconstructionDistribution dist = new ExponentialReconstructionDistribution("tanh");
            double negLogProb = dist.negLogProbability(x, distributionParams, average);
            INDArray exampleNegLogProb = dist.exampleNegLogProbability(x, distributionParams);
            assertArrayEquals(new int[] { minibatch, 1 }, exampleNegLogProb.shape());
            //Calculate the same thing, but using Apache Commons math
            double logProbSum = 0.0;
            for (int i = 0; i < minibatch; i++) {
                double exampleSum = 0.0;
                for (int j = 0; j < inputSize; j++) {
                    double gamma = gammas.getDouble(i, j);
                    double lambda = Math.exp(gamma);
                    double mean = 1.0 / lambda;
                    //Commons math uses mean = 1/lambda
                    ExponentialDistribution exp = new ExponentialDistribution(mean);
                    double xVal = x.getDouble(i, j);
                    double thisLogProb = exp.logDensity(xVal);
                    logProbSum += thisLogProb;
                    exampleSum += thisLogProb;
                }
                assertEquals(-exampleNegLogProb.getDouble(i), exampleSum, 1e-6);
            }
            double expNegLogProb;
            if (average) {
                expNegLogProb = -logProbSum / minibatch;
            } else {
                expNegLogProb = -logProbSum;
            }
            //                System.out.println(x);
            //                System.out.println(expNegLogProb + "\t" + logProb + "\t" + (logProb / expNegLogProb));
            assertEquals(expNegLogProb, negLogProb, 1e-6);
            //Also: check random sampling...
            int count = minibatch * inputSize;
            INDArray arr = Nd4j.linspace(-3, 3, count).reshape(minibatch, inputSize);
            INDArray sampleMean = dist.generateAtMean(arr);
            INDArray sampleRandom = dist.generateRandom(arr);
            for (int i = 0; i < minibatch; i++) {
                for (int j = 0; j < inputSize; j++) {
                    double d1 = sampleMean.getDouble(i, j);
                    double d2 = sampleRandom.getDouble(i, j);
                    assertTrue(d1 >= 0.0);
                    assertTrue(d2 >= 0.0);
                }
            }
        }
    }
}

Example 3 with ExponentialDistribution

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

the class BaseFunctionSolverTest method computeSCMOSWeights.

private static void computeSCMOSWeights(double[] weights, double[] noise) {
    // Per observation read noise.
    weights = new double[size * size];
    RandomGenerator randomGenerator = new Well19937c(42);
    ExponentialDistribution ed = new ExponentialDistribution(randomGenerator, variance, ExponentialDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
    for (int i = 0; i < weights.length; i++) {
        double pixelVariance = ed.sample();
        double pixelGain = Math.max(0.1, gain + randomGenerator.nextGaussian() * gainSD);
        // weights = var / g^2
        weights[i] = pixelVariance / (pixelGain * pixelGain);
    }
    // Convert to standard deviation for simulation
    noise = new double[weights.length];
    for (int i = 0; i < weights.length; i++) noise[i] = Math.sqrt(weights[i]);
}

Example 4 with ExponentialDistribution

use of org.apache.commons.math3.distribution.ExponentialDistribution in project incubator-systemml by apache.

the class ParameterizedBuiltin method computeFromDistribution.

/**
	 * Helper function to compute distribution-specific cdf (both lowertail and uppertail) and inverse cdf.
	 * 
	 * @param dcode probablility distribution code
	 * @param params map of parameters
	 * @param inverse true if inverse
	 * @return cdf or inverse cdf
	 * @throws MathArithmeticException if MathArithmeticException occurs
	 * @throws DMLRuntimeException if DMLRuntimeException occurs
	 */
private double computeFromDistribution(ProbabilityDistributionCode dcode, HashMap<String, String> params, boolean inverse) throws MathArithmeticException, DMLRuntimeException {
    // given value is "quantile" when inverse=false, and it is "probability" when inverse=true
    double val = Double.parseDouble(params.get("target"));
    boolean lowertail = true;
    if (params.get("lower.tail") != null) {
        lowertail = Boolean.parseBoolean(params.get("lower.tail"));
    }
    AbstractRealDistribution distFunction = null;
    switch(dcode) {
        case NORMAL:
            // default values for mean and sd
            double mean = 0.0, sd = 1.0;
            String mean_s = params.get("mean"), sd_s = params.get("sd");
            if (mean_s != null)
                mean = Double.parseDouble(mean_s);
            if (sd_s != null)
                sd = Double.parseDouble(sd_s);
            if (sd <= 0)
                throw new DMLRuntimeException("Standard deviation for Normal distribution must be positive (" + sd + ")");
            distFunction = new NormalDistribution(mean, sd);
            break;
        case EXP:
            // default value for 1/mean or rate
            double exp_rate = 1.0;
            if (params.get("rate") != null)
                exp_rate = Double.parseDouble(params.get("rate"));
            if (exp_rate <= 0) {
                throw new DMLRuntimeException("Rate for Exponential distribution must be positive (" + exp_rate + ")");
            }
            // For exponential distribution: mean = 1/rate
            distFunction = new ExponentialDistribution(1.0 / exp_rate);
            break;
        case CHISQ:
            if (params.get("df") == null) {
                throw new DMLRuntimeException("" + "Degrees of freedom must be specified for chi-squared distribution " + "(e.g., q=qchisq(0.5, df=20); p=pchisq(target=q, df=1.2))");
            }
            int df = UtilFunctions.parseToInt(params.get("df"));
            if (df <= 0) {
                throw new DMLRuntimeException("Degrees of Freedom for chi-squared distribution must be positive (" + df + ")");
            }
            distFunction = new ChiSquaredDistribution(df);
            break;
        case F:
            if (params.get("df1") == null || params.get("df2") == null) {
                throw new DMLRuntimeException("" + "Degrees of freedom must be specified for F distribution " + "(e.g., q = qf(target=0.5, df1=20, df2=30); p=pf(target=q, df1=20, df2=30))");
            }
            int df1 = UtilFunctions.parseToInt(params.get("df1"));
            int df2 = UtilFunctions.parseToInt(params.get("df2"));
            if (df1 <= 0 || df2 <= 0) {
                throw new DMLRuntimeException("Degrees of Freedom for F distribution must be positive (" + df1 + "," + df2 + ")");
            }
            distFunction = new FDistribution(df1, df2);
            break;
        case T:
            if (params.get("df") == null) {
                throw new DMLRuntimeException("" + "Degrees of freedom is needed to compute probabilities from t distribution " + "(e.g., q = qt(target=0.5, df=10); p = pt(target=q, df=10))");
            }
            int t_df = UtilFunctions.parseToInt(params.get("df"));
            if (t_df <= 0) {
                throw new DMLRuntimeException("Degrees of Freedom for t distribution must be positive (" + t_df + ")");
            }
            distFunction = new TDistribution(t_df);
            break;
        default:
            throw new DMLRuntimeException("Invalid distribution code: " + dcode);
    }
    double ret = Double.NaN;
    if (inverse) {
        // inverse cdf
        ret = distFunction.inverseCumulativeProbability(val);
    } else if (lowertail) {
        // cdf (lowertail)
        ret = distFunction.cumulativeProbability(val);
    } else {
        // cdf (upper tail)
        // TODO: more accurate distribution-specific computation of upper tail probabilities 
        ret = 1.0 - distFunction.cumulativeProbability(val);
    }
    return ret;
}