List of usage examples for org.apache.commons.math3.random RandomDataGenerator nextGaussian
public double nextGaussian(double mu, double sigma) throws NotStrictlyPositiveException
From source file:com.cloudera.oryx.ml.param.ContinuousAround.java
/** * @param rdg random number generator to use * @return a hyperparameter value chosen from Normal(around, step) *//*from w w w .j a va2 s .c om*/ @Override public Double getRandomValue(RandomDataGenerator rdg) { return rdg.nextGaussian(around, step); }
From source file:com.cloudera.oryx.ml.param.DiscreteAround.java
/** * @param rdg random number generator to use * @return a hyperparameter value chosen from Normal(around, step) and rounded to the nearest integer *//*from w w w. ja v a 2 s . c o m*/ @Override public Integer getRandomValue(RandomDataGenerator rdg) { return (int) Math.round(rdg.nextGaussian(around, step)); }
From source file:com.siemens.industrialbenchmark.util.RandomNumberExpectedValueTest.java
@Test public void testExpectedValues() { long seed = 0; Random rand = new Random(seed); RandomDataGenerator randomData = new RandomDataGenerator(); double uniformAverage = 0.0; double binomialAverage = 0.0; double normalAverage = 0.0; double exponentialAverage = 0.0; for (int i = 0; i < 1e6; i++) { // set current seed randomData.reSeed(seed);/* w w w. ja v a 2 s. c om*/ // draw random numbers double n = randomData.nextGaussian(0, 1); double u = randomData.nextUniform(0, 1); double b = randomData.nextBinomial(1, 0.5); double e = randomData.nextExponential(0.25); // average mean random number uniformAverage += (1. / (1. + i)) * (u - uniformAverage); binomialAverage += (1. / (1. + i)) * (b - binomialAverage); normalAverage += (1. / (1. + i)) * (n - normalAverage); exponentialAverage += (1. / (1. + i)) * (e - exponentialAverage); // draw new seed from global random generator seed = rand.nextLong(); } assertEquals(0.5, uniformAverage, 0.001); assertEquals(0.5, binomialAverage, 0.001); assertEquals(0.0, normalAverage, 0.001); assertEquals(0.25, exponentialAverage, 0.001); }
From source file:edu.cmu.tetrad.sem.LargeSemSimulator.java
public BoxDataSet simulateDataAcyclicConcurrent(int sampleSize) { int numVars = variableNodes.size(); setupModel(numVars);/* w w w.j ava2 s . c o m*/ System.out.println("Tier ordering"); final int[][] _parents = parents; final double[][] _coefs = coefs; // final double[][] _data = new double[sampleSize][numVars]; final double[][] _data = new double[numVars][sampleSize]; System.out.println("Starting simulation task"); // This random number generator is not thread safe, so we make a new one each time. RandomGenerator apacheGen = new Well19937a(new Date().getTime()); final RandomDataGenerator generator = new RandomDataGenerator(new SynchronizedRandomGenerator(apacheGen)); //Do the simulation. class SimulationTask extends RecursiveTask<Boolean> { private int chunk; private int from; private int to; public SimulationTask(int chunk, int from, int to) { this.chunk = chunk; this.from = from; this.to = to; } @Override protected Boolean compute() { RandomGenerator apacheGen = new Well44497b(generator.nextLong(0, Long.MAX_VALUE)); RandomDataGenerator generatorLocal = new RandomDataGenerator(apacheGen); if (to - from <= chunk) { for (int row = from; row < to; row++) { if (row % 100 == 0) out.println("Row " + row); for (int i = 0; i < tierIndices.length; i++) { int col = tierIndices[i]; double value = generatorLocal.nextGaussian(0, sqrt(errorVars[col])); for (int j = 0; j < _parents[col].length; j++) { int parent = _parents[col][j]; final double coef = _coefs[col][j]; final double v = _data[parent][row]; value += v * coef; if (Double.isNaN(value)) { throw new IllegalArgumentException(); } } value += means[col]; _data[col][row] = value; } } return true; } else { List<SimulationTask> simulationTasks = new ArrayList<SimulationTask>(); int mid = (to - from) / 2; simulationTasks.add(new SimulationTask(chunk, from, from + mid)); simulationTasks.add(new SimulationTask(chunk, from + mid, to)); invokeAll(simulationTasks); return true; } } } int chunk = 25; pool.invoke(new SimulationTask(chunk, 0, sampleSize)); return new BoxDataSet(new VerticalDoubleDataBox(_data), variableNodes); // return ColtDataSet.makeContinuousData(variableNodes, _data); }
From source file:gdsc.smlm.ij.plugins.pcpalm.PCPALMMolecules.java
private void runSimulation(boolean resultsAvailable) { if (resultsAvailable && !showSimulationDialog()) return;//ww w . jav a2 s . co m startLog(); log("Simulation parameters"); if (blinkingDistribution == 3) { log(" - Clusters = %d", nMolecules); log(" - Simulation size = %s um", Utils.rounded(simulationSize, 4)); log(" - Molecules/cluster = %s", Utils.rounded(blinkingRate, 4)); log(" - Blinking distribution = %s", BLINKING_DISTRIBUTION[blinkingDistribution]); log(" - p-Value = %s", Utils.rounded(p, 4)); } else { log(" - Molecules = %d", nMolecules); log(" - Simulation size = %s um", Utils.rounded(simulationSize, 4)); log(" - Blinking rate = %s", Utils.rounded(blinkingRate, 4)); log(" - Blinking distribution = %s", BLINKING_DISTRIBUTION[blinkingDistribution]); } log(" - Average precision = %s nm", Utils.rounded(sigmaS, 4)); log(" - Clusters simulation = " + CLUSTER_SIMULATION[clusterSimulation]); if (clusterSimulation > 0) { log(" - Cluster number = %s +/- %s", Utils.rounded(clusterNumber, 4), Utils.rounded(clusterNumberSD, 4)); log(" - Cluster radius = %s nm", Utils.rounded(clusterRadius, 4)); } final double nmPerPixel = 100; double width = simulationSize * 1000.0; // Allow a border of 3 x sigma for +/- precision //if (blinkingRate > 1) width -= 3 * sigmaS; RandomGenerator randomGenerator = new Well19937c( System.currentTimeMillis() + System.identityHashCode(this)); RandomDataGenerator dataGenerator = new RandomDataGenerator(randomGenerator); UniformDistribution dist = new UniformDistribution(null, new double[] { width, width, 0 }, randomGenerator.nextInt()); molecules = new ArrayList<Molecule>(nMolecules); // Create some dummy results since the calibration is required for later analysis results = new MemoryPeakResults(); results.setCalibration(new gdsc.smlm.results.Calibration(nmPerPixel, 1, 100)); results.setSource(new NullSource("Molecule Simulation")); results.begin(); int count = 0; // Generate a sequence of coordinates ArrayList<double[]> xyz = new ArrayList<double[]>((int) (nMolecules * 1.1)); Statistics statsRadius = new Statistics(); Statistics statsSize = new Statistics(); String maskTitle = TITLE + " Cluster Mask"; ByteProcessor bp = null; double maskScale = 0; // TODO - Add a fluctuations model to this. if (clusterSimulation > 0) { // Simulate clusters. // Note: In the Veatch et al. paper (Plos 1, e31457) correlation functions are built using circles // with small radii of 4-8 Arbitrary Units (AU) or large radii of 10-30 AU. A fluctuations model is // created at T = 1.075 Tc. It is not clear exactly how the particles are distributed. // It may be that a mask is created first using the model. The particles are placed on the mask using // a specified density. This simulation produces a figure to show either a damped cosine function // (circles) or an exponential (fluctuations). The number of particles in each circle may be randomly // determined just by density. The figure does not discuss the derivation of the cluster size // statistic. // // If this plugin simulation is run with a uniform distribution and blinking rate of 1 then the damped // cosine function is reproduced. The curve crosses g(r)=1 at a value equivalent to the average // distance to the centre-of-mass of each drawn cluster, not the input cluster radius parameter (which // is a hard upper limit on the distance to centre). final int maskSize = lowResolutionImageSize; int[] mask = null; maskScale = width / maskSize; // scale is in nm/pixel ArrayList<double[]> clusterCentres = new ArrayList<double[]>(); int totalSteps = 1 + (int) Math.ceil(nMolecules / clusterNumber); if (clusterSimulation == 2 || clusterSimulation == 3) { // Clusters are non-overlapping circles // Ensure the circles do not overlap by using an exclusion mask that accumulates // out-of-bounds pixels by drawing the last cluster (plus some border) on an image. When no // more pixels are available then stop generating molecules. // This is done by cumulatively filling a mask and using the MaskDistribution to select // a new point. This may be slow but it works. // TODO - Allow clusters of different sizes... mask = new int[maskSize * maskSize]; Arrays.fill(mask, 255); MaskDistribution maskDistribution = new MaskDistribution(mask, maskSize, maskSize, 0, maskScale, maskScale, randomGenerator); double[] centre; IJ.showStatus("Computing clusters mask"); int roiRadius = (int) Math.round((clusterRadius * 2) / maskScale); if (clusterSimulation == 3) { // Generate a mask of circles then sample from that. // If we want to fill the mask completely then adjust the total steps to be the number of // circles that can fit inside the mask. totalSteps = (int) (maskSize * maskSize / (Math.PI * Math.pow(clusterRadius / maskScale, 2))); } while ((centre = maskDistribution.next()) != null && clusterCentres.size() < totalSteps) { IJ.showProgress(clusterCentres.size(), totalSteps); // The mask returns the coordinates with the centre of the image at 0,0 centre[0] += width / 2; centre[1] += width / 2; clusterCentres.add(centre); // Fill in the mask around the centre to exclude any more circles that could overlap double cx = centre[0] / maskScale; double cy = centre[1] / maskScale; fillMask(mask, maskSize, (int) cx, (int) cy, roiRadius, 0); //log("[%.1f,%.1f] @ [%.1f,%.1f]", centre[0], centre[1], cx, cy); //Utils.display("Mask", new ColorProcessor(maskSize, maskSize, mask)); try { maskDistribution = new MaskDistribution(mask, maskSize, maskSize, 0, maskScale, maskScale, randomGenerator); } catch (IllegalArgumentException e) { // This can happen when there are no more non-zero pixels log("WARNING: No more room for clusters on the mask area (created %d of estimated %d)", clusterCentres.size(), totalSteps); break; } } IJ.showProgress(1); IJ.showStatus(""); } else { // Clusters are overlapping circles // Pick centres randomly from the distribution while (clusterCentres.size() < totalSteps) clusterCentres.add(dist.next()); } if (showClusterMask || clusterSimulation == 3) { // Show the mask for the clusters if (mask == null) mask = new int[maskSize * maskSize]; else Arrays.fill(mask, 0); int roiRadius = (int) Math.round((clusterRadius) / maskScale); for (double[] c : clusterCentres) { double cx = c[0] / maskScale; double cy = c[1] / maskScale; fillMask(mask, maskSize, (int) cx, (int) cy, roiRadius, 1); } if (clusterSimulation == 3) { // We have the mask. Now pick points at random from the mask. MaskDistribution maskDistribution = new MaskDistribution(mask, maskSize, maskSize, 0, maskScale, maskScale, randomGenerator); // Allocate each molecule position to a parent circle so defining clusters. int[][] clusters = new int[clusterCentres.size()][]; int[] clusterSize = new int[clusters.length]; for (int i = 0; i < nMolecules; i++) { double[] centre = maskDistribution.next(); // The mask returns the coordinates with the centre of the image at 0,0 centre[0] += width / 2; centre[1] += width / 2; xyz.add(centre); // Output statistics on cluster size and number. // TODO - Finding the closest cluster could be done better than an all-vs-all comparison double max = distance2(centre, clusterCentres.get(0)); int cluster = 0; for (int j = 1; j < clusterCentres.size(); j++) { double d2 = distance2(centre, clusterCentres.get(j)); if (d2 < max) { max = d2; cluster = j; } } // Assign point i to cluster centre[2] = cluster; if (clusterSize[cluster] == 0) { clusters[cluster] = new int[10]; } if (clusters[cluster].length <= clusterSize[cluster]) { clusters[cluster] = Arrays.copyOf(clusters[cluster], (int) (clusters[cluster].length * 1.5)); } clusters[cluster][clusterSize[cluster]++] = i; } // Generate real cluster size statistics for (int j = 0; j < clusterSize.length; j++) { final int size = clusterSize[j]; if (size == 0) continue; statsSize.add(size); if (size == 1) { statsRadius.add(0); continue; } // Find centre of cluster and add the distance to each point double[] com = new double[2]; for (int n = 0; n < size; n++) { double[] xy = xyz.get(clusters[j][n]); for (int k = 0; k < 2; k++) com[k] += xy[k]; } for (int k = 0; k < 2; k++) com[k] /= size; for (int n = 0; n < size; n++) { double dx = xyz.get(clusters[j][n])[0] - com[0]; double dy = xyz.get(clusters[j][n])[1] - com[1]; statsRadius.add(Math.sqrt(dx * dx + dy * dy)); } } } if (showClusterMask) { bp = new ByteProcessor(maskSize, maskSize); for (int i = 0; i < mask.length; i++) if (mask[i] != 0) bp.set(i, 128); Utils.display(maskTitle, bp); } } // Use the simulated cluster centres to create clusters of the desired size if (clusterSimulation == 1 || clusterSimulation == 2) { for (double[] clusterCentre : clusterCentres) { int clusterN = (int) Math.round( (clusterNumberSD > 0) ? dataGenerator.nextGaussian(clusterNumber, clusterNumberSD) : clusterNumber); if (clusterN < 1) continue; //double[] clusterCentre = dist.next(); if (clusterN == 1) { // No need for a cluster around a point xyz.add(clusterCentre); statsRadius.add(0); statsSize.add(1); } else { // Generate N random points within a circle of the chosen cluster radius. // Locate the centre-of-mass and the average distance to the centre. double[] com = new double[3]; int j = 0; while (j < clusterN) { // Generate a random point within a circle uniformly // http://stackoverflow.com/questions/5837572/generate-a-random-point-within-a-circle-uniformly double t = 2.0 * Math.PI * randomGenerator.nextDouble(); double u = randomGenerator.nextDouble() + randomGenerator.nextDouble(); double r = clusterRadius * ((u > 1) ? 2 - u : u); double x = r * Math.cos(t); double y = r * Math.sin(t); double[] xy = new double[] { clusterCentre[0] + x, clusterCentre[1] + y }; xyz.add(xy); for (int k = 0; k < 2; k++) com[k] += xy[k]; j++; } // Add the distance of the points from the centre of the cluster. // Note this does not account for the movement due to precision. statsSize.add(j); if (j == 1) { statsRadius.add(0); } else { for (int k = 0; k < 2; k++) com[k] /= j; while (j > 0) { double dx = xyz.get(xyz.size() - j)[0] - com[0]; double dy = xyz.get(xyz.size() - j)[1] - com[1]; statsRadius.add(Math.sqrt(dx * dx + dy * dy)); j--; } } } } } } else { // Random distribution for (int i = 0; i < nMolecules; i++) xyz.add(dist.next()); } // The Gaussian sigma should be applied so the overall distance from the centre // ( sqrt(x^2+y^2) ) has a standard deviation of sigmaS? final double sigma1D = sigmaS / Math.sqrt(2); // Show optional histograms StoredDataStatistics intraDistances = null; StoredDataStatistics blinks = null; if (showHistograms) { int capacity = (int) (xyz.size() * blinkingRate); intraDistances = new StoredDataStatistics(capacity); blinks = new StoredDataStatistics(capacity); } Statistics statsSigma = new Statistics(); for (int i = 0; i < xyz.size(); i++) { int nOccurrences = getBlinks(dataGenerator, blinkingRate); if (showHistograms) blinks.add(nOccurrences); final int size = molecules.size(); // Get coordinates in nm final double[] moleculeXyz = xyz.get(i); if (bp != null && nOccurrences > 0) { bp.putPixel((int) Math.round(moleculeXyz[0] / maskScale), (int) Math.round(moleculeXyz[1] / maskScale), 255); } while (nOccurrences-- > 0) { final double[] localisationXy = Arrays.copyOf(moleculeXyz, 2); // Add random precision if (sigma1D > 0) { final double dx = dataGenerator.nextGaussian(0, sigma1D); final double dy = dataGenerator.nextGaussian(0, sigma1D); localisationXy[0] += dx; localisationXy[1] += dy; if (!dist.isWithinXY(localisationXy)) continue; // Calculate mean-squared displacement statsSigma.add(dx * dx + dy * dy); } final double x = localisationXy[0]; final double y = localisationXy[1]; molecules.add(new Molecule(x, y, i, 1)); // Store in pixels float[] params = new float[7]; params[Gaussian2DFunction.X_POSITION] = (float) (x / nmPerPixel); params[Gaussian2DFunction.Y_POSITION] = (float) (y / nmPerPixel); results.add(i + 1, (int) x, (int) y, 0, 0, 0, params, null); } if (molecules.size() > size) { count++; if (showHistograms) { int newCount = molecules.size() - size; if (newCount == 1) { // No intra-molecule distances //intraDistances.add(0); continue; } // Get the distance matrix between these molecules double[][] matrix = new double[newCount][newCount]; for (int ii = size, x = 0; ii < molecules.size(); ii++, x++) { for (int jj = size + 1, y = 1; jj < molecules.size(); jj++, y++) { final double d2 = molecules.get(ii).distance2(molecules.get(jj)); matrix[x][y] = matrix[y][x] = d2; } } // Get the maximum distance for particle linkage clustering of this molecule double max = 0; for (int x = 0; x < newCount; x++) { // Compare to all-other molecules and get the minimum distance // needed to join at least one double linkDistance = Double.POSITIVE_INFINITY; for (int y = 0; y < newCount; y++) { if (x == y) continue; if (matrix[x][y] < linkDistance) linkDistance = matrix[x][y]; } // Check if this is larger if (max < linkDistance) max = linkDistance; } intraDistances.add(Math.sqrt(max)); } } } results.end(); if (bp != null) Utils.display(maskTitle, bp); // Used for debugging //System.out.printf(" * Molecules = %d (%d activated)\n", xyz.size(), count); //if (clusterSimulation > 0) // System.out.printf(" * Cluster number = %s +/- %s. Radius = %s +/- %s\n", // Utils.rounded(statsSize.getMean(), 4), Utils.rounded(statsSize.getStandardDeviation(), 4), // Utils.rounded(statsRadius.getMean(), 4), Utils.rounded(statsRadius.getStandardDeviation(), 4)); log("Simulation results"); log(" * Molecules = %d (%d activated)", xyz.size(), count); log(" * Blinking rate = %s", Utils.rounded((double) molecules.size() / xyz.size(), 4)); log(" * Precision (Mean-displacement) = %s nm", (statsSigma.getN() > 0) ? Utils.rounded(Math.sqrt(statsSigma.getMean()), 4) : "0"); if (showHistograms) { if (intraDistances.getN() == 0) { log(" * Mean Intra-Molecule particle linkage distance = 0 nm"); log(" * Fraction of inter-molecule particle linkage @ 0 nm = 0 %%"); } else { plot(blinks, "Blinks/Molecule", true); double[][] intraHist = plot(intraDistances, "Intra-molecule particle linkage distance", false); // Determine 95th and 99th percentile int p99 = intraHist[0].length - 1; double limit1 = 0.99 * intraHist[1][p99]; double limit2 = 0.95 * intraHist[1][p99]; while (intraHist[1][p99] > limit1 && p99 > 0) p99--; int p95 = p99; while (intraHist[1][p95] > limit2 && p95 > 0) p95--; log(" * Mean Intra-Molecule particle linkage distance = %s nm (95%% = %s, 99%% = %s, 100%% = %s)", Utils.rounded(intraDistances.getMean(), 4), Utils.rounded(intraHist[0][p95], 4), Utils.rounded(intraHist[0][p99], 4), Utils.rounded(intraHist[0][intraHist[0].length - 1], 4)); if (distanceAnalysis) { performDistanceAnalysis(intraHist, p99); } } } if (clusterSimulation > 0) { log(" * Cluster number = %s +/- %s", Utils.rounded(statsSize.getMean(), 4), Utils.rounded(statsSize.getStandardDeviation(), 4)); log(" * Cluster radius = %s +/- %s nm (mean distance to centre-of-mass)", Utils.rounded(statsRadius.getMean(), 4), Utils.rounded(statsRadius.getStandardDeviation(), 4)); } }