List of usage examples for org.apache.commons.math3.distribution RealDistribution getNumericalMean
double getNumericalMean();
From source file:com.wormsim.utils.Utils.java
/** * Returns a string representation of the provided distribution. TODO: Make * this complete TODO: Make this compatible with custom distributions (or just * more complex ones).//w w w . java2s .c o m * * @param dist The distribution to translate * * @return The distribution as a string. */ public static String realDistributionToString(RealDistribution dist) { if (dist instanceof ConstantRealDistribution) { return Double.toString(dist.getNumericalMean()); } else if (dist instanceof UniformRealDistribution) { return "Uniform(" + dist.getSupportLowerBound() + "," + dist.getSupportUpperBound() + ")"; } else if (dist instanceof NormalDistribution) { NormalDistribution dist2 = (NormalDistribution) dist; return "Normal(" + dist2.getMean() + "," + dist2.getStandardDeviation() + ")"; } else { return dist.toString(); } }
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>/* w w w . ja v a 2 s . co 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:org.nmdp.ngs.reads.GeneratePairedEndReads.java
/** * Generate paired-end next generation sequencing (NGS/HTS) reads. * * @param reference reference, must not be null * @param variant FASTQ variant, must not be null * @param random random, must not be null * @param length length distribution, must not be null * @param insertSize insert size distribution, must not be null * @param quality quality strategy, must not be null * @param coverage coverage strategy, must not be null * @param mutationRate mutation rate, must be between <code>0.0</code> and <code>1.0</code>, inclusive * @param mutation mutation strategy, must not be null * @param first first appendable, must not be null * @param second second appendable, must not be null * @param writer FASTQ writer, must not be null *//*from w w w .ja v a2s. c om*/ public GeneratePairedEndReads(final Sequence reference, final FastqVariant variant, final RandomGenerator random, final RealDistribution length, final RealDistribution insertSize, final QualityStrategy quality, final CoverageStrategy coverage, final double mutationRate, final MutationStrategy mutation, final Appendable first, final Appendable second, final FastqWriter writer) { checkNotNull(reference, "reference must not be null"); checkNotNull(variant, "variant must not be null"); checkNotNull(random, "random must not be null"); checkNotNull(length, "length must not be null"); checkNotNull(insertSize, "insertSize must not be null"); checkNotNull(quality, "quality must not be null"); checkNotNull(coverage, "coverage must not be null"); checkArgument(mutationRate >= 0.0d, "mutationRate must be greater than or equal to 0.0d"); checkArgument(mutationRate <= 1.0d, "mutationRate must be less than or equal to 1.0d"); checkNotNull(mutation, "mutation must not be null"); checkNotNull(first, "first must not be null"); checkNotNull(second, "second must not be null"); checkNotNull(writer, "writer must not be null"); this.reference = reference; this.variant = variant; this.random = random; this.length = length; this.insertSize = insertSize; this.quality = quality; this.coverage = coverage; this.mutationRate = mutationRate; this.mutation = mutation; this.first = first; this.second = second; this.writer = writer; int flanking = (int) (length.getNumericalMean() + length.getNumericalVariance() + insertSize.getNumericalMean() + insertSize.getNumericalVariance()); location = new UniformIntegerDistribution(this.random, 1 - flanking, this.reference.length() + flanking); }
From source file:org.nmdp.ngs.reads.GenerateReads.java
/** * Generate next generation sequencing (NGS/HTS) reads. * * @param reference reference, must not be null * @param variant FASTQ variant, must not be null * @param random random, must not be null * @param length length distribution, must not be null * @param quality quality strategy, must not be null * @param coverage coverage strategy, must not be null * @param mutationRate mutation rate, must be between <code>0.0</code> and <code>1.0</code>, inclusive * @param mutation mutation strategy, must not be null * @param appendable appendable, must not be null * @param writer FASTQ writer, must not be null *//* w w w . j a va2 s .com*/ public GenerateReads(final Sequence reference, final FastqVariant variant, final RandomGenerator random, final RealDistribution length, final QualityStrategy quality, final CoverageStrategy coverage, final double mutationRate, final MutationStrategy mutation, final Appendable appendable, final FastqWriter writer) { checkNotNull(reference, "reference must not be null"); checkNotNull(variant, "variant must not be null"); checkNotNull(random, "random must not be null"); checkNotNull(length, "length must not be null"); checkNotNull(quality, "quality must not be null"); checkNotNull(coverage, "coverage must not be null"); checkArgument(mutationRate >= 0.0d, "mutationRate must be greater than or equal to 0.0d"); checkArgument(mutationRate <= 1.0d, "mutationRate must be less than or equal to 1.0d"); checkNotNull(mutation, "mutation must not be null"); checkNotNull(appendable, "appendable must not be null"); checkNotNull(writer, "writer must not be null"); this.reference = reference; this.variant = variant; this.random = random; this.length = length; this.quality = quality; this.coverage = coverage; this.mutationRate = mutationRate; this.mutation = mutation; this.appendable = appendable; this.writer = writer; int flanking = (int) (length.getNumericalMean() + length.getNumericalVariance()); location = new UniformIntegerDistribution(this.random, 1 - flanking, this.reference.length() + flanking); }
From source file:org.nmdp.ngs.tools.GenerateBed.java
/** * Generate records in BED format./*from ww w.ja va 2 s. com*/ * * @param bedFile output BED file, if any * @param n number of BED records to generate, must be at least zero * @param size chromosome size, must be at least zero * @param chrom chromosome name, must not be null * @param random random generator, must not be null * @param length length distribution, must not be null */ public GenerateBed(final File bedFile, final int n, final int size, final String chrom, final RandomGenerator random, final RealDistribution length) { checkArgument(n >= 0, "n must be at least zero"); checkArgument(size >= 0, "size must be at least zero"); checkNotNull(chrom); checkNotNull(random); checkNotNull(length); this.bedFile = bedFile; this.n = n; this.size = size; this.chrom = chrom; this.random = random; this.length = length; int flanking = (int) (length.getNumericalMean() + length.getNumericalVariance()); location = new UniformIntegerDistribution(this.random, 1 - flanking, size + flanking); }