Java Code Examples for org.apache.commons.math3.distribution.PoissonDistribution#sample()

The following examples show how to use org.apache.commons.math3.distribution.PoissonDistribution#sample() . You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example 1
Source File: RandSPInstruction.java    From systemds with Apache License 2.0 5 votes vote down vote up 
@Override
public Iterator<Double> call(SampleTask t)
		throws Exception {

	long st = t.range_start;
	long end = Math.min(t.range_start+_partitionSize, _maxValue);
	ArrayList<Double> retList = new ArrayList<>();
	
	if ( _frac == 1.0 ) 
	{
		for(long i=st; i <= end; i++) 
			retList.add((double)i);
	}
	else 
	{
		if(_replace) 
		{
			PoissonDistribution pdist = new PoissonDistribution( (_frac > 0.0 ? _frac :1.0) );
			for(long i=st; i <= end; i++)
			{
				int count = pdist.sample();
				while(count > 0) {
					retList.add((double)i);
					count--;
				}
			}
		}
		else 
		{
			Random rnd = new Random(t.seed);
			for(long i=st; i <=end; i++) 
				if ( rnd.nextDouble() < _frac )
					retList.add((double) i);
		}
	}
	return retList.iterator();
}
Example 2
Source File: RandSPInstruction.java    From systemds with Apache License 2.0 5 votes vote down vote up 
@Override
public Iterator<Double> call(SampleTask t)
		throws Exception {

	long st = t.range_start;
	long end = Math.min(t.range_start+_partitionSize, _maxValue);
	ArrayList<Double> retList = new ArrayList<>();
	
	if ( _frac == 1.0 ) 
	{
		for(long i=st; i <= end; i++) 
			retList.add((double)i);
	}
	else 
	{
		if(_replace) 
		{
			PoissonDistribution pdist = new PoissonDistribution( (_frac > 0.0 ? _frac :1.0) );
			for(long i=st; i <= end; i++)
			{
				int count = pdist.sample();
				while(count > 0) {
					retList.add((double)i);
					count--;
				}
			}
		}
		else 
		{
			Random rnd = new Random(t.seed);
			for(long i=st; i <=end; i++) 
				if ( rnd.nextDouble() < _frac )
					retList.add((double) i);
		}
	}
	return retList.iterator();
}
Example 3
Source File: RandomDataGeneratorTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up 
@Test
/**
 * MATH-720
 */
public void testReseed() {
    PoissonDistribution x = new PoissonDistribution(3.0);
    x.reseedRandomGenerator(0);
    final double u = x.sample();
    PoissonDistribution y = new PoissonDistribution(3.0);
    y.reseedRandomGenerator(0);
    Assert.assertEquals(u, y.sample(), 0);
}
Example 4
Source File: RandomDataTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up 
@Test
/**
 * MATH-720
 */
public void testReseed() {
    PoissonDistribution x = new PoissonDistribution(3.0);
    x.reseedRandomGenerator(0);
    final double u = x.sample();
    PoissonDistribution y = new PoissonDistribution(3.0);
    y.reseedRandomGenerator(0);
    Assert.assertEquals(u, y.sample(), 0);
}
Example 5
Source File: RandomDataTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up 
@Test
/**
 * MATH-720
 */
public void testReseed() {
    PoissonDistribution x = new PoissonDistribution(3.0);
    x.reseedRandomGenerator(0);
    final double u = x.sample();
    PoissonDistribution y = new PoissonDistribution(3.0);
    y.reseedRandomGenerator(0);
    Assert.assertEquals(u, y.sample(), 0);
}
Example 6
Source File: RandomDataGeneratorTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up 
@Test
/**
 * MATH-720
 */
public void testReseed() {
    PoissonDistribution x = new PoissonDistribution(3.0);
    x.reseedRandomGenerator(0);
    final double u = x.sample();
    PoissonDistribution y = new PoissonDistribution(3.0);
    y.reseedRandomGenerator(0);
    Assert.assertEquals(u, y.sample(), 0);
}
Example 7
Source File: RandomDataTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up 
@Test
/**
 * MATH-720
 */
public void testReseed() {
    PoissonDistribution x = new PoissonDistribution(3.0);
    x.reseedRandomGenerator(0);
    final double u = x.sample();
    PoissonDistribution y = new PoissonDistribution(3.0);
    y.reseedRandomGenerator(0);
    Assert.assertEquals(u, y.sample(), 0);
}
Example 8
Source File: RandomDataGeneratorTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up 
@Test
/**
 * MATH-720
 */
public void testReseed() {
    PoissonDistribution x = new PoissonDistribution(3.0);
    x.reseedRandomGenerator(0);
    final double u = x.sample();
    PoissonDistribution y = new PoissonDistribution(3.0);
    y.reseedRandomGenerator(0);
    Assert.assertEquals(u, y.sample(), 0);
}
Example 9
Source File: RandomDataGeneratorTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up 
@Test
/**
 * MATH-720
 */
public void testReseed() {
    PoissonDistribution x = new PoissonDistribution(3.0);
    x.reseedRandomGenerator(0);
    final double u = x.sample();
    PoissonDistribution y = new PoissonDistribution(3.0);
    y.reseedRandomGenerator(0);
    Assert.assertEquals(u, y.sample(), 0);
}
Example 10
Source File: RandomDataGeneratorTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up 
@Test
/**
 * MATH-720
 */
public void testReseed() {
    PoissonDistribution x = new PoissonDistribution(3.0);
    x.reseedRandomGenerator(0);
    final double u = x.sample();
    PoissonDistribution y = new PoissonDistribution(3.0);
    y.reseedRandomGenerator(0);
    Assert.assertEquals(u, y.sample(), 0);
}
Example 11
Source File: ChangeFinder2DTest.java    From incubator-hivemall with Apache License 2.0 4 votes vote down vote up 
public void testSota5D() throws HiveException {
    final int DIM = 5;
    final int EXAMPLES = 20001;

    final Double[] x = new Double[DIM];
    final List<Double> xList = Arrays.asList(x);

    Parameters params = new Parameters();
    params.set(LossFunction.logloss);
    params.r1 = 0.01d;
    params.k = 10;
    params.T1 = 10;
    params.T2 = 10;
    PrimitiveObjectInspector oi = PrimitiveObjectInspectorFactory.javaDoubleObjectInspector;
    ListObjectInspector listOI = ObjectInspectorFactory.getStandardListObjectInspector(oi);
    final ChangeFinder2D cf = new ChangeFinder2D(params, listOI);
    final double[] outScores = new double[2];

    RandomGenerator rng1 = new Well19937c(31L);
    final UniformIntegerDistribution uniform = new UniformIntegerDistribution(rng1, 0, 10);
    RandomGenerator rng2 = new Well19937c(41L);
    final PoissonDistribution poissonEvent = new PoissonDistribution(rng2, 1000.d,
        PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);
    final StringBuilder buf = new StringBuilder(256);

    println("# time x0 x1 x2 x3 x4 mean0 mean1 mean2 mean3 mean4 outlier change");
    FIN: for (int i = 0; i < EXAMPLES;) {
        int len = poissonEvent.sample();
        double data[][] = new double[DIM][len];
        double mean[] = new double[DIM];
        double sd[] = new double[DIM];
        for (int j = 0; j < DIM; j++) {
            mean[j] = uniform.sample() * 5.d;
            sd[j] = uniform.sample() / 10.d * 5.d + 1.d;
            if (i % 5 == 0) {
                mean[j] += 50.d;
            }
            NormalDistribution normDist =
                    new NormalDistribution(new Well19937c(i + j), mean[j], sd[j]);
            data[j] = normDist.sample(len);
            data[j][len / (j + 2) + DIM % (j + 1)] = mean[j] + (j + 4) * sd[j];
        }
        for (int j = 0; j < len; j++) {
            if (i >= EXAMPLES) {
                break FIN;
            }
            x[0] = data[0][j];
            x[1] = data[1][j];
            x[2] = data[2][j];
            x[3] = data[3][j];
            x[4] = data[4][j];
            cf.update(xList, outScores);
            buf.append(i)
               .append(' ')
               .append(x[0].doubleValue())
               .append(' ')
               .append(x[1].doubleValue())
               .append(' ')
               .append(x[2].doubleValue())
               .append(' ')
               .append(x[3].doubleValue())
               .append(' ')
               .append(x[4].doubleValue())
               .append(' ')
               .append(mean[0])
               .append(' ')
               .append(mean[1])
               .append(' ')
               .append(mean[2])
               .append(' ')
               .append(mean[3])
               .append(' ')
               .append(mean[4])
               .append(' ')
               .append(outScores[0])
               .append(' ')
               .append(outScores[1]);
            println(buf.toString());
            StringUtils.clear(buf);
            i++;
        }
    }
}
Example 12
Source File: AlleleFractionSimulatedData.java    From gatk-protected with BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up 
public AlleleFractionSimulatedData(final double averageHetsPerSegment, final int numSegments,
        final double averageDepth, final double biasMean, final double biasVariance, final double outlierProbability) {
    rng.setSeed(RANDOM_SEED);
    this.numSegments = numSegments;
    final AlleleFractionState.MinorFractions minorFractions = new AlleleFractionState.MinorFractions(numSegments);
    final List<AllelicCount> alleleCounts = new ArrayList<>();
    final List<SimpleInterval> segments = new ArrayList<>();

    final PoissonDistribution segmentLengthGenerator = makePoisson(rng, averageHetsPerSegment);
    final PoissonDistribution readDepthGenerator = makePoisson(rng, averageDepth);
    final UniformRealDistribution minorFractionGenerator = new UniformRealDistribution(rng, 0.0, 0.5);

    //translate to ApacheCommons' parametrization of the gamma distribution
    final double gammaShape = biasMean * biasMean / biasVariance;
    final double gammaScale = biasVariance / biasMean;
    final GammaDistribution biasGenerator = new GammaDistribution(rng, gammaShape, gammaScale);

    //put each segment on its own chromosome and sort by lexicographical order
    final List<String> chromosomes = IntStream.range(0, numSegments).mapToObj(Integer::toString).collect(Collectors.toList());
    Collections.sort(chromosomes);

    for (final String chromosome : chromosomes) {
        // calculate the range of het indices for this segment
        final int numHetsInSegment = Math.max(MIN_HETS_PER_SEGMENT, segmentLengthGenerator.sample());

        final double minorFraction = minorFractionGenerator.sample();
        minorFractions.add(minorFraction);

        //we will put all the hets in this segment/chromosome at loci 1, 2, 3 etc
        segments.add(new SimpleInterval(chromosome, 1, numHetsInSegment + 1));
        for (int het = 1; het < numHetsInSegment + 1; het++) {
            final double bias = biasGenerator.sample();

            //flip a coin to decide alt minor (alt fraction = minor fraction) or ref minor (alt fraction = 1 - minor fraction)
            final boolean isAltMinor = rng.nextDouble() < 0.5;
            final double altFraction =  isAltMinor ? minorFraction : 1 - minorFraction;

            //the probability of an alt read is the alt fraction modified by the bias or, in the case of an outlier, random
            final double pAlt;
            if (rng.nextDouble() < outlierProbability) {
                truePhases.add(AlleleFractionIndicator.OUTLIER);
                pAlt = rng.nextDouble();
            } else {
                truePhases.add(isAltMinor ? AlleleFractionIndicator.ALT_MINOR : AlleleFractionIndicator.REF_MINOR);
                pAlt = altFraction / (altFraction + (1 - altFraction) * bias);
            }

            final int numReads = readDepthGenerator.sample();
            final int numAltReads = new BinomialDistribution(rng, numReads, pAlt).sample();
            final int numRefReads = numReads - numAltReads;
            alleleCounts.add(new AllelicCount(new SimpleInterval(chromosome, het, het), numRefReads, numAltReads));
        }
    }

    final Genome genome = new Genome(TRIVIAL_TARGETS, alleleCounts);
    segmentedGenome = new SegmentedGenome(segments, genome);
    trueState = new AlleleFractionState(biasMean, biasVariance, outlierProbability, minorFractions);
}
Example 13
Source File: AlleleFractionSimulatedData.java    From gatk with BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up 
AlleleFractionSimulatedData(final SampleLocatableMetadata metadata,
                            final AlleleFractionGlobalParameters globalParameters,
                            final int numSegments,
                            final double averageHetsPerSegment,
                            final double averageDepth,
                            final RandomGenerator rng) {
    final AlleleFractionState.MinorFractions minorFractions = new AlleleFractionState.MinorFractions(numSegments);
    final List<AllelicCount> allelicCounts = new ArrayList<>();
    final List<SimpleInterval> segments = new ArrayList<>();

    final PoissonDistribution segmentLengthGenerator = makePoisson(rng, averageHetsPerSegment);
    final PoissonDistribution readDepthGenerator = makePoisson(rng, averageDepth);
    final UniformRealDistribution minorFractionGenerator = new UniformRealDistribution(rng, 0.0, 0.5);

    final double meanBias = globalParameters.getMeanBias();
    final double biasVariance = globalParameters.getBiasVariance();
    final double outlierProbability = globalParameters.getOutlierProbability();

    //translate to ApacheCommons' parametrization of the gamma distribution
    final double gammaShape = meanBias * meanBias / biasVariance;
    final double gammaScale = biasVariance / meanBias;
    final GammaDistribution biasGenerator = new GammaDistribution(rng, gammaShape, gammaScale);

    //put each segment on its own chromosome and sort in sequence-dictionary order
    final List<String> chromosomes = IntStream.range(0, numSegments)
            .mapToObj(i -> metadata.getSequenceDictionary().getSequence(i).getSequenceName())
            .collect(Collectors.toList());

    for (final String chromosome : chromosomes) {
        // calculate the range of het indices for this segment
        final int numHetsInSegment = Math.max(MIN_HETS_PER_SEGMENT, segmentLengthGenerator.sample());

        final double minorFraction = minorFractionGenerator.sample();
        minorFractions.add(minorFraction);

        //we will put all the hets in this segment/chromosome at loci 1, 2, 3 etc
        segments.add(new SimpleInterval(chromosome, 1, numHetsInSegment));
        for (int het = 1; het < numHetsInSegment + 1; het++) {
            final double bias = biasGenerator.sample();

            //flip a coin to decide alt minor (alt fraction = minor fraction) or ref minor (alt fraction = 1 - minor fraction)
            final boolean isAltMinor = rng.nextDouble() < 0.5;
            final double altFraction =  isAltMinor ? minorFraction : 1 - minorFraction;

            //the probability of an alt read is the alt fraction modified by the bias or, in the case of an outlier, random
            final double pAlt;
            if (rng.nextDouble() < outlierProbability) {
                pAlt = rng.nextDouble();
            } else {
                pAlt = altFraction / (altFraction + (1 - altFraction) * bias);
            }

            final int numReads = readDepthGenerator.sample();
            final int numAltReads = new BinomialDistribution(rng, numReads, pAlt).sample();
            final int numRefReads = numReads - numAltReads;
            allelicCounts.add(new AllelicCount(new SimpleInterval(chromosome, het, het), numRefReads, numAltReads));
        }
    }

    data = new AlleleFractionSegmentedData(
            new AllelicCountCollection(metadata, allelicCounts),
            new SimpleIntervalCollection(metadata, segments));
    trueState = new AlleleFractionState(meanBias, biasVariance, outlierProbability, minorFractions);
}