use of org.apache.commons.math3.stat.descriptive.moment.Mean in project gatk by broadinstitute.
the class SliceSamplerUnitTest method testInitialPointOutOfRange.
@Test(expectedExceptions = IllegalArgumentException.class)
public void testInitialPointOutOfRange() {
rng.setSeed(RANDOM_SEED);
final double mean = 5.;
final double standardDeviation = 0.75;
final NormalDistribution normalDistribution = new NormalDistribution(mean, standardDeviation);
final Function<Double, Double> normalLogPDF = normalDistribution::logDensity;
final double xInitial = -10.;
final double xMin = 0.;
final double xMax = 1.;
final double width = 0.5;
final SliceSampler normalSampler = new SliceSampler(rng, normalLogPDF, xMin, xMax, width);
normalSampler.sample(xInitial);
}
use of org.apache.commons.math3.stat.descriptive.moment.Mean in project gatk by broadinstitute.
the class SliceSamplerUnitTest method testSliceSamplingOfNormalDistribution.
/**
* Test slice sampling of a normal distribution. Checks that input mean and standard deviation are recovered
* by 10000 samples to a relative error of 0.5% and 2%, respectively.
*/
@Test
public void testSliceSamplingOfNormalDistribution() {
rng.setSeed(RANDOM_SEED);
final double mean = 5.;
final double standardDeviation = 0.75;
final NormalDistribution normalDistribution = new NormalDistribution(mean, standardDeviation);
final Function<Double, Double> normalLogPDF = normalDistribution::logDensity;
final double xInitial = 1.;
final double xMin = Double.NEGATIVE_INFINITY;
final double xMax = Double.POSITIVE_INFINITY;
final double width = 0.5;
final int numSamples = 10000;
final SliceSampler normalSampler = new SliceSampler(rng, normalLogPDF, xMin, xMax, width);
final double[] samples = Doubles.toArray(normalSampler.sample(xInitial, numSamples));
final double sampleMean = new Mean().evaluate(samples);
final double sampleStandardDeviation = new StandardDeviation().evaluate(samples);
Assert.assertEquals(relativeError(sampleMean, mean), 0., 0.005);
Assert.assertEquals(relativeError(sampleStandardDeviation, standardDeviation), 0., 0.02);
}
use of org.apache.commons.math3.stat.descriptive.moment.Mean in project gatk by broadinstitute.
the class AdaptiveMetropolisSamplerUnitTest method testBeta.
@Test
public void testBeta() {
final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(RANDOM_SEED));
for (final double a : Arrays.asList(10, 20, 30)) {
for (final double b : Arrays.asList(10, 20, 30)) {
final double theoreticalMean = a / (a + b);
final double theoreticalVariance = a * b / ((a + b) * (a + b) * (a + b + 1));
//Note: this is the theoretical standard deviation of the sample mean given uncorrelated
//samples. The sample mean will have a greater variance here because samples are correlated.
final double standardDeviationOfMean = Math.sqrt(theoreticalVariance / NUM_SAMPLES);
final Function<Double, Double> logPDF = x -> (a - 1) * Math.log(x) + (b - 1) * Math.log(1 - x);
final AdaptiveMetropolisSampler sampler = new AdaptiveMetropolisSampler(INITIAL_BETA_GUESS, INITIAL_STEP_SIZE, 0, 1);
final List<Double> samples = sampler.sample(rng, logPDF, NUM_SAMPLES, NUM_BURN_IN_STEPS);
final double sampleMean = samples.stream().mapToDouble(x -> x).average().getAsDouble();
final double sampleMeanSquare = samples.stream().mapToDouble(x -> x * x).average().getAsDouble();
final double sampleVariance = (sampleMeanSquare - sampleMean * sampleMean) * NUM_SAMPLES / (NUM_SAMPLES - 1);
Assert.assertEquals(sampleMean, theoreticalMean, 10 * standardDeviationOfMean);
Assert.assertEquals(sampleVariance, theoreticalVariance, 10e-4);
}
}
}
use of org.apache.commons.math3.stat.descriptive.moment.Mean in project gatk by broadinstitute.
the class AdaptiveMetropolisSamplerUnitTest method testGaussian.
@Test
public void testGaussian() {
final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(RANDOM_SEED));
for (final double theoreticalMean : Arrays.asList(0)) {
for (final double precision : Arrays.asList(1.0)) {
final double variance = 1 / precision;
//Note: this is the theoretical standard deviation of the sample mean given uncorrelated
//samples. The sample mean will have a greater variance here because samples are correlated.
final double standardDeviationOfMean = Math.sqrt(variance / NUM_SAMPLES);
final Function<Double, Double> logPDF = x -> -(precision / 2) * (x - theoreticalMean) * (x - theoreticalMean);
final AdaptiveMetropolisSampler sampler = new AdaptiveMetropolisSampler(INITIAL_GAUSSIAN_GUESS, INITIAL_STEP_SIZE, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);
final List<Double> samples = sampler.sample(rng, logPDF, NUM_SAMPLES, NUM_BURN_IN_STEPS);
final double sampleMean = samples.stream().mapToDouble(x -> x).average().getAsDouble();
final double sampleMeanSquare = samples.stream().mapToDouble(x -> x * x).average().getAsDouble();
final double sampleVariance = (sampleMeanSquare - sampleMean * sampleMean) * NUM_SAMPLES / (NUM_SAMPLES - 1);
Assert.assertEquals(sampleMean, theoreticalMean, 6 * standardDeviationOfMean);
Assert.assertEquals(sampleVariance, variance, variance / 10);
}
}
}
use of org.apache.commons.math3.stat.descriptive.moment.Mean in project gatk-protected by broadinstitute.
the class CNLOHCaller method calculateVarianceOfCopyNeutralSegmentMeans.
/**
* Attempt to get an idea of segment mean variance near copy neutral.
*
* @param segments Never {@code null}
* @return variance of segment mean (in CR space) of segments that are "close enough" to copy neutral.
* Zero if no segments are "close enough"
*/
private double calculateVarianceOfCopyNeutralSegmentMeans(final List<ACNVModeledSegment> segments, final double meanBiasInCR) {
Utils.nonNull(segments);
// Only consider values "close enough" to copy neutral (CR == 1).
final double neutralCR = 1 + meanBiasInCR;
final double[] neutralSegmentMeans = segments.stream().mapToDouble(ACNVModeledSegment::getSegmentMeanInCRSpace).filter(m -> Math.abs(m - neutralCR) < CLOSE_ENOUGH_TO_COPY_NEUTRAL_IN_CR).toArray();
return new Variance().evaluate(neutralSegmentMeans);
}
Aggregations