org.apache.commons.math3.distribution.PoissonDistribution Java Examples

/**
 * Calculates k, based on Zero-truncated Poisson distribution.
 * https://en.wikipedia.org/wiki/Zero-truncated_Poisson_distribution
 * 
 * @param lambda
 * @return
 */
private int calculateKZeroPoisson(double lambda) {
    
    final double threshold = 1d - this.significanceLevel;
    final PoissonDistribution distribution = new PoissonDistribution(lambda);
    final double v2 = 1d - distribution.probability(0);
    int counter = 1;
    double value = 0d;
    while (value < threshold) {
        // value2 += ((Math.pow(lambda, counter)) / (Math.exp(lambda) - 1)) * ArithmeticUtils.factorial(counter);
        value += distribution.probability(counter) / v2;
        counter++;
        // Break if the estimated k is equal or greater than the given k, as this makes no sense.
        if (counter >= this.k) {
            // We are 100% sure that the dataset fulfills k-map
            value = 1d;
            break;
        }
    }
    this.type1Error = 1d - value;
    return counter;
}
 

/**
 * Calculates k, based on Poisson distribution.
 * @param lambda
 * @return
 */
private int calculateKPoisson(double lambda) {
    
    final double threshold = 1d - this.significanceLevel;
    final PoissonDistribution distribution = new PoissonDistribution(lambda);
    int counter = 0;
    double value = 0;
    while (value < threshold) {
        // value += (Math.pow(lambda, counter) * Math.exp(-lambda)) / ArithmeticUtils.factorial(counter);
        value = distribution.cumulativeProbability(counter);
        counter++;
        // Break if the estimated k is equal or greater than the given k, as this makes no sense.
        if (counter >= this.k) {
            // We are 100% sure that the dataset fulfills k-map
            value = 1d;
            break;
        }
    }
    this.type1Error = 1d - value;
    return counter + 1;
}
 

/**
 * Composes a isd descriptor adding random spaces b
 * @param baseName the distribution name.
 * @param mean the mean for the distribution.
 * @param stddev the stddev for the distribution.
 * @param spacesDistr poisson distribution of random spaces to add in different parts of the descriptor.
 * @return never {@code null}.
 */
private String composeDescriptionString(final String baseName, final double mean, final double stddev, PoissonDistribution spacesDistr) {
    return StringUtils.repeat(' ', spacesDistr.sample()) +
            baseName +
            StringUtils.repeat(' ', spacesDistr.sample()) +
            "(" +
            StringUtils.repeat(' ', spacesDistr.sample()) +
            mean +
            StringUtils.repeat(' ', spacesDistr.sample()) +
            ',' +
            StringUtils.repeat(' ', spacesDistr.sample()) +
            stddev +
            StringUtils.repeat(' ', spacesDistr.sample()) +
            ')' +
            StringUtils.repeat(' ', spacesDistr.sample());
}
 

/**
 * Test the sampler functions at the upper bound on the mean.
 */
@Test
public void testSamplerAtUpperBounds() {
    final double mean = SUPPORTED_UPPER_BOUND;

    // Test some ranges for the cumulative probability
    final PoissonDistribution pd = new PoissonDistribution(null, mean,
        PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);

    final FixedRNG rng = new FixedRNG(0);
    final SharedStateDiscreteSampler sampler = KempSmallMeanPoissonSampler.of(rng, mean);

    // Lower bound should be zero
    testSample(rng, sampler, 0, 0, 0);

    // Upper bound should exceed 99.99% of the range
    testSample(rng, sampler, 1, pd.inverseCumulativeProbability(0.9999), Integer.MAX_VALUE);

    // A sample from within the cumulative probability should be within the expected range
    for (int i = 1; i < 10; i++) {
        final double p = i * 0.1;
        final int lower = pd.inverseCumulativeProbability(p - 0.01);
        final int upper = pd.inverseCumulativeProbability(p + 0.01);
        testSample(rng, sampler, p, lower, upper);
    }
}
 

static double probabilityDriverVariant(long firstVariantTypeCount, long secondVariantTypeCount,
        @NotNull final DndsDriverImpactLikelihood firstLikelihood, @NotNull final DndsDriverImpactLikelihood secondLikelihood) {
    if (firstLikelihood.equals(secondLikelihood)) {
        return probabilityDriverVariantSameImpact(1, firstVariantTypeCount, firstLikelihood);
    }

    double lambda1 = firstVariantTypeCount * firstLikelihood.pVariantNonDriverFactor();
    double lambda2 = secondVariantTypeCount * secondLikelihood.pVariantNonDriverFactor();
    if (Doubles.isZero(lambda1) || Doubles.isZero(lambda2)) {
        return Math.max(probabilityDriverVariant(firstVariantTypeCount, firstLikelihood),
                probabilityDriverVariant(secondVariantTypeCount, secondLikelihood));
    }

    final double pDriver = Math.max(firstLikelihood.pDriver(), secondLikelihood.pDriver());
    final double pVariantNonDriver1 = 1 - new PoissonDistribution(lambda1).cumulativeProbability(0);
    final double pVariantNonDriver2 = 1 - new PoissonDistribution(lambda2).cumulativeProbability(0);
    final double pVariantNonDriver = pVariantNonDriver1 * pVariantNonDriver2;

    return pDriver / (pDriver + pVariantNonDriver * (1 - pDriver));
}
 

/** {@inheritDoc} */
@Override public BootstrappedDatasetPartition build(
    LearningEnvironment env,
    Iterator<UpstreamEntry<K, V>> upstreamData,
    long upstreamDataSize,
    EmptyContext ctx) {

    BootstrappedVector[] dataset = new BootstrappedVector[Math.toIntExact(upstreamDataSize)];

    int cntr = 0;

    PoissonDistribution poissonDistribution = new PoissonDistribution(
        new Well19937c(env.randomNumbersGenerator().nextLong()),
        subsampleSize,
        PoissonDistribution.DEFAULT_EPSILON,
        PoissonDistribution.DEFAULT_MAX_ITERATIONS);

    while (upstreamData.hasNext()) {
        UpstreamEntry<K, V> nextRow = upstreamData.next();
        LabeledVector<Double> vecAndLb = preprocessor.apply(nextRow.getKey(), nextRow.getValue());
        Vector features = vecAndLb.features();
        Double lb = vecAndLb.label();
        int[] repetitionCounters = new int[samplesCnt];
        Arrays.setAll(repetitionCounters, i -> poissonDistribution.sample());
        dataset[cntr++] = new BootstrappedVector(features, lb, repetitionCounters);
    }

    return new BootstrappedDatasetPartition(dataset);
}
 

@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();
}
 

/** {@inheritDoc} */
@Override public Stream<UpstreamEntry> transform(Stream<UpstreamEntry> upstream) {
    PoissonDistribution poisson = new PoissonDistribution(
        new Well19937c(seed),
        subsampleRatio,
        PoissonDistribution.DEFAULT_EPSILON,
        PoissonDistribution.DEFAULT_MAX_ITERATIONS);

    return upstream.sequential().flatMap(en -> Stream.generate(() -> en).limit(poisson.sample()));
}
 

/**
 * Create a poisson sampler which can sample elements with replacement.
 *
 * @param fraction The expected count of each element.
 * @param seed     Random number generator seed for internal PoissonDistribution.
 */
public PoissonSampler(double fraction, long seed) {
	Preconditions.checkArgument(fraction >= 0, "fraction should be positive.");
	this.fraction = fraction;
	if (this.fraction > 0) {
		this.poissonDistribution = new PoissonDistribution(fraction);
		this.poissonDistribution.reseedRandomGenerator(seed);
	}
	this.random = new XORShiftRandom(seed);
}
 

@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);
}
 

@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);
}
 

@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);
}
 

@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);
}
 

@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);
}
 

@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);
}
 

@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);
}
 

@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);
}
 

/**
 * Create a poisson sampler which can sample elements with replacement.
 *
 * @param fraction The expected count of each element.
 * @param seed     Random number generator seed for internal PoissonDistribution.
 */
public PoissonSampler(double fraction, long seed) {
	Preconditions.checkArgument(fraction >= 0, "fraction should be positive.");
	this.fraction = fraction;
	if (this.fraction > 0) {
		this.poissonDistribution = new PoissonDistribution(fraction);
		this.poissonDistribution.reseedRandomGenerator(seed);
	}
	this.random = new XORShiftRandom(seed);
}
 

/**
 * Create a poisson sampler which can sample elements with replacement.
 *
 * @param fraction The expected count of each element.
 */
public PoissonSampler(double fraction) {
	Preconditions.checkArgument(fraction >= 0, "fraction should be non-negative.");
	this.fraction = fraction;
	if (this.fraction > 0) {
		this.poissonDistribution = new PoissonDistribution(fraction);
	}
	this.random = new XORShiftRandom();
}
 

@Inject
PoissonPoisonedContextStartHandler(
    @Parameter(CrashProbability.class) final double lambda, final Clock clock) {

  this.clock = clock;
  this.timeToCrash = new PoissonDistribution(lambda * 1000).sample();

  LOG.log(Level.INFO,
      "Created Poisson poison injector with prescribed dose: {0}. Crash in {1} msec.",
      new Object[]{lambda, this.timeToCrash});
}
 

@Inject
PoissonPoisonedTaskStartHandler(
    @Parameter(CrashProbability.class) final double lambda, final Clock clock) {

  this.clock = clock;
  this.timeToCrash = new PoissonDistribution(lambda * 1000).sample();

  LOG.log(Level.INFO,
      "Created Poisson poison injector with prescribed dose: {0}. Crash in {1} msec.",
      new Object[]{lambda, this.timeToCrash});
}
 

@Override
public Object doWork(Object first) throws IOException{
  if(null == first){
    throw new IOException(String.format(Locale.ROOT,"Invalid expression %s - null found for the first value",toExpression(constructingFactory)));
  }

  Number mean = (Number)first;

  return new PoissonDistribution(mean.doubleValue());
}
 

/**
 * Create a poisson sampler which can sample elements with replacement.
 *
 * @param fraction The expected count of each element.
 */
public PoissonSampler(double fraction) {
	Preconditions.checkArgument(fraction >= 0, "fraction should be non-negative.");
	this.fraction = fraction;
	if (this.fraction > 0) {
		this.poissonDistribution = new PoissonDistribution(fraction);
	}
	this.random = new XORShiftRandom();
}
 

/**
 * Create a poisson sampler which can sample elements with replacement.
 *
 * @param fraction The expected count of each element.
 * @param seed     Random number generator seed for internal PoissonDistribution.
 */
public PoissonSampler(double fraction, long seed) {
	Preconditions.checkArgument(fraction >= 0, "fraction should be positive.");
	this.fraction = fraction;
	if (this.fraction > 0) {
		this.poissonDistribution = new PoissonDistribution(fraction);
		this.poissonDistribution.reseedRandomGenerator(seed);
	}
	this.random = new XORShiftRandom(seed);
}
 

/**
 * Create a poisson sampler which can sample elements with replacement.
 *
 * @param fraction The expected count of each element.
 */
public PoissonSampler(double fraction) {
	Preconditions.checkArgument(fraction >= 0, "fraction should be non-negative.");
	this.fraction = fraction;
	if (this.fraction > 0) {
		this.poissonDistribution = new PoissonDistribution(fraction);
	}
	this.random = new XORShiftRandom();
}
 

@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();
}
 

@Override
public Double map(Row row) throws Exception {
    double label = (Double) row.getField(numFeature);
    double weight = (Double) row.getField(numFeature + 1);
    double pred = (Double) row.getField(numFeature + 3);

    PoissonDistribution distribution = new PoissonDistribution(pred);
    return weight * Math.log(distribution.probability((int) label));
}
 

/**
 * Test sampling from a Poisson distribution.
 */
@Test
public void testPoissonSamples() {
    final double mean = 3.14;
    final PoissonDistribution dist = new PoissonDistribution(null, mean,
        PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);
    final int maxN = dist.inverseCumulativeProbability(1 - 1e-6);
    final double[] expected = new double[maxN];
    for (int i = 0; i < expected.length; i++) {
        expected[i] = dist.probability(i);
    }
    checkSamples(expected, 1.0);
}
 

/**
 * Test sampling from a Poisson distribution.
 */
@Test
public void testPoissonSamples() {
    final double mean = 3.14;
    final PoissonDistribution dist = new PoissonDistribution(null, mean,
        PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);
    final int maxN = dist.inverseCumulativeProbability(1 - 1e-6);
    double[] expected = new double[maxN];
    for (int i = 0; i < expected.length; i++) {
        expected[i] = dist.probability(i);
    }
    checkSamples(expected);
}
 

/**
 * Calculates the likelihood of different sex genotypes for a given sample index
 *
 * @param sampleIndex sample index
 * @return an instance of {@link SexGenotypeData}
 */
private SexGenotypeData calculateSexGenotypeData(final int sampleIndex) {
    final int[] allosomalReadCounts = Arrays.stream(processedReadCounts.getColumnOnSpecifiedTargets(sampleIndex,
            allosomalTargetList, true)).mapToInt(n -> (int)n).toArray();
    final double readDepth = getSampleReadDepthFromAutosomalTargets(sampleIndex);

    logger.debug("Read depth for " + processedReadCounts.columnNames().get(sampleIndex) + ": " + readDepth);
    final List<Double> logLikelihoods = new ArrayList<>();
    final List<String> sexGenotypesList = new ArrayList<>();
    final int numAllosomalTargets = allosomalTargetList.size();

    for (final String genotypeName : sexGenotypeIdentifiers) {
        /* calculate log likelihood */
        final int[] currentAllosomalTargetPloidies = allosomalTargetPloidies.get(genotypeName);
        final double[] poissonParameters = IntStream.range(0, numAllosomalTargets)
                .mapToDouble(i -> readDepth *
                        (currentAllosomalTargetPloidies[i] > 0 ? currentAllosomalTargetPloidies[i]
                                : baselineMappingErrorProbability))
                .toArray();
        final double currentLogLikelihood = IntStream.range(0, numAllosomalTargets)
                .mapToDouble(i -> {
                    final PoissonDistribution pois = new PoissonDistribution(poissonParameters[i]);
                    return pois.logProbability(allosomalReadCounts[i]);
                }).sum();
        sexGenotypesList.add(genotypeName);
        logLikelihoods.add(currentLogLikelihood);
    }

    /* infer the most likely sex genotype */
    final Integer[] indices = new Integer[sexGenotypesList.size()];
    IntStream.range(0, sexGenotypesList.size()).forEach(i -> indices[i] = i);
    Arrays.sort(indices, (li, ri) -> Double.compare(logLikelihoods.get(ri), logLikelihoods.get(li)));

    return new SexGenotypeData(processedReadCounts.columnNames().get(sampleIndex),
            sexGenotypesList.get(indices[0]), sexGenotypesList,
            logLikelihoods.stream().mapToDouble(d -> d).toArray());
}