List of usage examples for org.apache.commons.math3.distribution GammaDistribution GammaDistribution
public GammaDistribution(RandomGenerator rng, double shape, double scale, double inverseCumAccuracy) throws NotStrictlyPositiveException
From source file:gdsc.smlm.model.ImageModel.java
/** * Simulate an image of fluorophores. The total amount of time a fluorophore is on (i.e. sum of tOn) is used to * create the signal strength using the specified correlation. * <p>/*from w ww . j a va 2 s .c o m*/ * A second set of fluorophores, independent of the first, are generated using * {@link #createFluorophores(List, int)}. The correlated on-times will be created by combining the times using the * provided correlation (r): * * <pre> * // X1 : Fluorophore total on-times * // X2 : Independently generated fluorophore total on-times * // r : correlation * a = sqrt(1 - r * r) * newX = (r * X1 + a * X2) / (r + a) * </pre> * * The signal is proportional to newly generated on-times (newX) with an average of the specified photon budget. * <p> * The photon budget can either be distributed evenly over the fluorophore lifetime or per frame (see * {@link #isPhotonBudgetPerFrame()}). Each frame signal output will be subject to Poisson sampling. * <p> * If the input correlation is zero then the number of photons will be sampled from the configured photon * distribution or, if this is null, from a gamma distribution with a mean of the specified photon budget and the * shape parameter defined by {@link #setPhotonShapeParameter(double)}. * <p> * A random fraction of the fluorophores will move using a random walk with the diffusion coefficient defined in the * compound. * * @param compoundFluorophores * The compounds containing the fluorophores * @param fixedFraction * The fraction of molecules that are fixed * @param maxFrames * The maximum frame for the simulation * @param photonBudget * the average number of photons per frame/lifetime; see {@link #isPhotonBudgetPerFrame()} * @param correlation * The correlation between the number of photons and the total on-time of the fluorophore * @param rotate * Rotate the molecule per frame * @return the localisations */ public List<LocalisationModel> createImage(List<CompoundMoleculeModel> compoundFluorophores, double fixedFraction, int maxFrames, double photonBudget, double correlation, boolean rotate) { List<LocalisationModel> localisations = new ArrayList<LocalisationModel>(); // Extract the fluorophores in all the compounds ArrayList<FluorophoreSequenceModel> fluorophores = new ArrayList<FluorophoreSequenceModel>( compoundFluorophores.size()); for (CompoundMoleculeModel c : compoundFluorophores) { for (int i = c.getSize(); i-- > 0;) if (c.getMolecule(i) instanceof FluorophoreSequenceModel) { fluorophores.add((FluorophoreSequenceModel) c.getMolecule(i)); } } final int nMolecules = fluorophores.size(); // Check the correlation bounds. // Correlation for photons per frame verses total on time should be negative. // Correlation for total photons verses total on time should be positive. double r; if (photonBudgetPerFrame) r = (correlation < -1) ? -1 : (correlation > 0) ? 0 : correlation; else r = (correlation < 0) ? 0 : (correlation > 1) ? 1 : correlation; // Create a photon budget for each fluorophore double[] photons = new double[nMolecules]; // Generate a second set of on times using the desired correlation if (r != 0) { // Observations on real data show: // - there is a weak positive correlation between total photons and time // - There is a weak negative correlation between photons per frame and total on-time // How to generate a random correlation: // http://www.uvm.edu/~dhowell/StatPages/More_Stuff/CorrGen.html // http://stats.stackexchange.com/questions/13382/how-to-define-a-distribution-such-that-draws-from-it-correlate-with-a-draw-from // Used for debugging the correlation //StoredDataStatistics onTime = new StoredDataStatistics(); //StoredDataStatistics allPhotons = new StoredDataStatistics(); //PearsonsCorrelation c = new PearsonsCorrelation(); // Create a second set of fluorophores. This is used to generate the correlated photon data List<FluorophoreSequenceModel> fluorophores2 = new ArrayList<FluorophoreSequenceModel>(); while (fluorophores2.size() < fluorophores.size()) { FluorophoreSequenceModel f = createFluorophore(0, new double[] { 0, 0, 0 }, maxFrames); if (f != null) fluorophores2.add(f); } final double a = Math.sqrt(1 - r * r); // Q - How to generate a negative correlation? // Generate a positive correlation then invert the data and shift to the same range boolean invert = (r < 0); r = Math.abs(r); StoredDataStatistics correlatedOnTime = new StoredDataStatistics(); for (int i = 0; i < nMolecules; i++) { final double X = getTotalOnTime(fluorophores.get(i)); final double Z = getTotalOnTime(fluorophores2.get(i)); //onTime.add(X); correlatedOnTime.add((r * X + a * Z) / (r + a)); } if (invert) { final double min = correlatedOnTime.getStatistics().getMin(); final double range = correlatedOnTime.getStatistics().getMax() - min; StoredDataStatistics newCorrelatedOnTime = new StoredDataStatistics(); for (double d : correlatedOnTime.getValues()) { newCorrelatedOnTime.add(range - d + min); } correlatedOnTime = newCorrelatedOnTime; } // Get the average on time from the correlated sample. // Using the population value (tOn * (1+nBlinks)) over-estimates the on time. final double averageTotalTOn = correlatedOnTime.getMean(); // Now create the localisations double[] correlatedOnTimes = correlatedOnTime.getValues(); for (int i = 0; i < nMolecules; i++) { // Generate photons using the correlated on time data double p = photonBudget * correlatedOnTimes[i] / averageTotalTOn; //final double X = getTotalOnTime(fluorophores.get(i)); //double p = (X > 0) ? photonBudget * correlatedOnTimes[i] / X : 0; //double p = (X > 0) ? randomGenerator.nextGamma(photonBudget, correlatedOnTimes[i] / X) : 0; //double p = randomGenerator.nextGamma(photonBudget, correlatedOnTimes[i] / X); //allPhotons.add(p); photons[i] = p; } //System.out.printf("t = %f, p = %f, R = %f\n", averageTotalTOn, allPhotons.getMean(), // c.correlation(onTime.getValues(), allPhotons.getValues())); } else { // Sample from the provided distribution. Do not over-write the class level distribution to allow // running again with a different shape parameter / photon budget. RealDistribution photonDistribution = getPhotonDistribution(); final double photonScale; if (photonDistribution == null) { photonScale = 1; // Model using a gamma distribution if we have a shape parameter if (photonShapeParameter > 0) { final double scaleParameter = photonBudget / photonShapeParameter; photonDistribution = new GammaDistribution(random, photonShapeParameter, scaleParameter, ExponentialDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY); } } else { // Ensure the custom distribution is scaled to the correct photon budget photonScale = photonBudget / photonDistribution.getNumericalMean(); } if (photonDistribution == null) { // No distribution so use the same number for all Arrays.fill(photons, photonBudget); } else { // Generate photons sampling from the photon budget for (int i = 0; i < nMolecules; i++) { photons[i] = photonDistribution.sample() * photonScale; } } } int photonIndex = 0; for (CompoundMoleculeModel c : compoundFluorophores) { final boolean fixed = (random.nextDouble() < fixedFraction); photonIndex += createLocalisation(c, localisations, fixed, maxFrames, photons, photonIndex, !fixed && rotate); } sortByTime(localisations); return localisations; }
From source file:msi.gaml.operators.Stats.java
@operator(value = "gamma_rnd", can_be_const = false, type = IType.LIST, category = { IOperatorCategory.STATISTICAL }, concept = { IConcept.STATISTIC, IConcept.CLUSTERING }) @doc(value = "returns a random value from a gamma distribution with specified values of the shape and scale parameters", examples = { @example(value = "gamma_rnd(10.0,5.0)", isExecutable = false) }) public static Double OpGammaDist(final IScope scope, final Double shape, final Double scale) throws GamaRuntimeException { final GammaDistribution dist = new GammaDistribution(scope.getRandom().getGenerator(), shape, scale, GammaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY); return dist.sample(); }
From source file:org.jreserve.jrlib.util.random.RndGammaTest.java
@Before public void setUp() { rnd = new RndGamma(new JavaRandom(SEED)); double scale = VARIANCE / MEAN; double shape = MEAN / scale; RandomGenerator rg = new JDKRandomGenerator(); rg.setSeed(SEED);/*from ww w .j av a 2s .c om*/ gd = new GammaDistribution(rg, shape, scale, GammaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY); }