import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
public void normallyRandomInitialParams(double mu, double sigma) {
	RandomGenerator rg = RandomGeneratorFactory.createRandomGenerator(new Random(System.currentTimeMillis()));
	_W = new double[layerNum - 1][][];
	for (int l = 0; l < layerNum - 1; l++) {
		_W[l] = new double[activationLayerNum[l + 1]][];
		for (int i = 0; i < activationLayerNum[l + 1]; i++) {
			_W[l][i] = new double[activationLayerNum[l]];
			for (int j = 0; j < activationLayerNum[l]; j++) {
				// 将原始为mu为0.5，3sigma为0.5的正态分布，转移到mu，sigma为此处设置的正态分布的随机值
				_W[l][i][j] = (rg.nextGaussian() - 0.5 + mu) * 6.0 * sigma;
			}
		}
	}

	_B = new double[layerNum - 1];

	for (int l = 0; l < layerNum - 1; l++) {
		_B[l] = (rg.nextGaussian() - 0.5 + mu) * 6.0 * sigma;
	}
}
 

import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
@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);
        }
    }
}
 

import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
@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);
        }
    }
}
 

import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
@Test
public void testChromosomesOnDifferentSegments() {
    final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(563));

    final double[] trueMinorAlleleFractions = new double[] {0.12, 0.32, 0.5};
    final double trueMemoryLength = 1e5;
    final AlleleFractionGlobalParameters trueParams = new AlleleFractionGlobalParameters(1.0, 0.01, 0.01);

    // randomly set positions
    final int chainLength = 100;
    final List<SimpleInterval> positions = CopyRatioSegmenterUnitTest.randomPositions("chr1", chainLength, rng, trueMemoryLength/4);
    positions.addAll(CopyRatioSegmenterUnitTest.randomPositions("chr2", chainLength, rng, trueMemoryLength/4));
    positions.addAll(CopyRatioSegmenterUnitTest.randomPositions("chr3", chainLength, rng, trueMemoryLength/4));

    final int trueState = 2;    //fix everything to the same state 2
    final List<Double> minorAlleleFractionSequence = Collections.nCopies(positions.size(), trueMinorAlleleFractions[trueState]);

    final AllelicCountCollection counts = generateCounts(minorAlleleFractionSequence, positions, rng, trueParams);

    final AlleleFractionSegmenter segmenter = new AlleleFractionSegmenter(10, counts, AllelicPanelOfNormals.EMPTY_PON);
    final List<ModeledSegment> segments = segmenter.getModeledSegments();

    //check that each chromosome has at least one segment
    final int numDifferentContigsInSegments = (int) segments.stream().map(ModeledSegment::getContig).distinct().count();
    Assert.assertEquals(numDifferentContigsInSegments, 3);
}
 

import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
private Object[][] getUnivariateGaussianTargets(final double sigma) {
    Random rng = new Random(337);
    final RandomGenerator randomGenerator = RandomGeneratorFactory.createRandomGenerator(rng);
    NormalDistribution n = new NormalDistribution(randomGenerator, 1, sigma);
    final int numDataPoints = 10000;
    final int numEventPoints = 2000;

    final List<ReadCountRecord.SingleSampleRecord> targetList = new ArrayList<>();
    for (int i = 0; i < (numDataPoints - numEventPoints); i++){
        targetList.add(new ReadCountRecord.SingleSampleRecord(new Target("arbitrary_name", new SimpleInterval("chr1", 100 + 2*i, 101 + 2 * i)), n.sample()));
    }
    for (int i = (numDataPoints - numEventPoints); i < numDataPoints; i++){
        targetList.add(new ReadCountRecord.SingleSampleRecord(new Target("arbitrary_name", new SimpleInterval("chr1", 100 + 2 * i, 101 + 2 * i)), 0.5 + n.sample()));
    }

    HashedListTargetCollection<ReadCountRecord.SingleSampleRecord> targets = new HashedListTargetCollection<>(targetList);

    List<ModeledSegment> segments = new ArrayList<>();
    segments.add(new ModeledSegment(new SimpleInterval("chr1", 100, 16050), 8000, 1));
    segments.add(new ModeledSegment(new SimpleInterval("chr1", 16100, 20200), 2000, 1.5));

    return new Object [] []{ {targets, segments}};
}
 

import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
@Test
public void testMCMC() {
    final double variance = 0.01;
    final double outlierProbability = 0.05;
    final int numSegments = 100;
    final double averageIntervalsPerSegment = 100.;
    final int numSamples = 150;
    final int numBurnIn = 50;
    final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(RANDOM_SEED));

    final SampleLocatableMetadata metadata = new SimpleSampleLocatableMetadata(
            "test-sample",
            new SAMSequenceDictionary(IntStream.range(0, numSegments)
                    .mapToObj(i -> new SAMSequenceRecord("chr" + i + 1, 10000))
                    .collect(Collectors.toList())));
    final CopyRatioSimulatedData simulatedData = new CopyRatioSimulatedData(
            metadata, variance, outlierProbability, numSegments, averageIntervalsPerSegment, rng);

    final CopyRatioModeller modeller = new CopyRatioModeller(simulatedData.getData().getCopyRatios(), simulatedData.getData().getSegments());
    modeller.fitMCMC(numSamples, numBurnIn);

    assertCopyRatioPosteriorCenters(modeller, simulatedData);
}
 

import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
public BM(int numClusters, int dimension, long randomSeed) {
    this.numClusters = numClusters;
    this.dimension = dimension;
    this.distributions = new BernoulliDistribution[numClusters][dimension];
    this.mixtureCoefficients = new double[numClusters];
    Arrays.fill(mixtureCoefficients,1.0/numClusters);
    this.logMixtureCoefficients = new double[numClusters];
    Arrays.fill(logMixtureCoefficients,Math.log(1.0/numClusters));
    Random random = new Random(randomSeed);
    RandomGenerator randomGenerator = RandomGeneratorFactory.createRandomGenerator(random);
    UniformRealDistribution uniform = new UniformRealDistribution(randomGenerator, 0.25,0.75);
    for (int k=0;k<numClusters;k++){
        for (int d=0;d<dimension;d++){
            double p = uniform.sample();
            distributions[k][d] = new BernoulliDistribution(p);
        }
    }
    this.logClusterConditioinalForEmpty = new double[numClusters];
    updateLogClusterConditioinalForEmpty();

    this.names = new ArrayList<>(dimension);
    for (int d=0;d<dimension;d++){
        names.add(""+d);
    }
}
 

import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
default List<S> generateHiddenStateChain(final List<T> positions) {
    final RandomGenerator rg = RandomGeneratorFactory.createRandomGenerator(new Random(RANDOM_SEED_FOR_CHAIN_GENERATION));
    final List<S> hiddenStates = hiddenStates();
    final List<S> result = new ArrayList<>(positions.size());

    final S initialState = GATKProtectedMathUtils.randomSelect(hiddenStates, s -> Math.exp(logPriorProbability(s, positions.get(0))), rg);
    result.add(initialState);

    IntStream.range(1, positions.size()).forEach(n ->
        result.add(GATKProtectedMathUtils.randomSelect(hiddenStates,
                s -> Math.exp(logTransitionProbability(result.get(n-1), positions.get(n - 1), s, positions.get(n))), rg)));
    return result;
}
 

import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
/**
 * Creates a new model instance.
 *
 * @param eventStartProbability the probability per base pair of a transition from neutral to a CNV
 * @param meanEventSize the expectation of the distance between consecutive targets in an event
 * @param deletionMeanShift the deletion depth of coverage negative shift.
 * @param duplicationMeaShift the duplication depth of coverage positive shift.
 * @throws IllegalArgumentException if any of the model parameters has an invalid value.
 */
public XHMMModel(final double eventStartProbability,
                 final double meanEventSize,
                 final double deletionMeanShift,
                 final double duplicationMeaShift) {
    super(new XHMMEmissionProbabilityCalculator(
            deletionMeanShift,
            duplicationMeaShift,
            EMISSION_SD, RandomGeneratorFactory.createRandomGenerator(new Random(RANDOM_SEED))),
            eventStartProbability,
            meanEventSize);
}
 

import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
/** <p>Produces a Gaussian mixture model based on the difference between targets and segment means.</p>
 * <p>Filters targets to populations where more than the minProportion lie in a single segment.</p>
 * <p>Returns null if no pass filtering.  Please note that in these cases,
 * in the rest of this class, we use this to assume that a GMM is not a good model.</p>
 *
 * @param segments  -- segments with segment mean in log2 copy ratio space
 * @param targets -- targets with a log2 copy ratio estimate
 * @param minProportion -- minimum proportion of all targets that a given segment must have in order to be used
 *                      in the evaluation
 * @param numComponents -- number of components to use in the GMM.  Usually, this is 2.
 * @return  never {@code null}.  Fitting result with indications whether it converged or was even attempted.
 */
private MixtureMultivariateNormalFitResult retrieveGaussianMixtureModelForFilteredTargets(final List<ModeledSegment> segments,
                                                                                          final TargetCollection<ReadCountRecord.SingleSampleRecord> targets, double minProportion, int numComponents){

    // For each target in a segment that contains enough targets, normalize the difference against the segment mean
    //  and collapse the filtered targets into the copy ratio estimates.
    final List<Double> filteredTargetsSegDiff = getNumProbeFilteredTargetList(segments, targets, minProportion);

    if (filteredTargetsSegDiff.size() < numComponents) {
        return new MixtureMultivariateNormalFitResult(null, false, false);
    }

    // Assume that Apache Commons wants data points in the first dimension.
    // Note that second dimension of length 2 (instead of 1) is to wrok around funny Apache commons API.
    final double[][] filteredTargetsSegDiff2d = new double[filteredTargetsSegDiff.size()][2];

    // Convert the filtered targets into 2d array (even if second dimension is length 1).  The second dimension is
    //  uncorrelated Gaussian.  This is only to get around funny API in Apache Commons, which will throw an
    //  exception if the length of the second dimension is < 2
    final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(RANDOM_SEED));
    final NormalDistribution nd = new NormalDistribution(rng, 0, .1);
    for (int i = 0; i < filteredTargetsSegDiff.size(); i++) {
        filteredTargetsSegDiff2d[i][0] = filteredTargetsSegDiff.get(i);
        filteredTargetsSegDiff2d[i][1] = nd.sample();
    }

    final MixtureMultivariateNormalDistribution estimateEM0 = MultivariateNormalMixtureExpectationMaximization.estimate(filteredTargetsSegDiff2d, numComponents);
    final MultivariateNormalMixtureExpectationMaximization multivariateNormalMixtureExpectationMaximization = new MultivariateNormalMixtureExpectationMaximization(filteredTargetsSegDiff2d);

    try {
        multivariateNormalMixtureExpectationMaximization.fit(estimateEM0);
    } catch (final MaxCountExceededException | ConvergenceException | SingularMatrixException e) {
        // We are done, we cannot make a fitting.  We should return a result that we attempted a fitting, but it
        //  did not converge.  Include the model as it was when the exception was thrown.
        return new MixtureMultivariateNormalFitResult(multivariateNormalMixtureExpectationMaximization.getFittedModel(), false, true);
    }
    return new MixtureMultivariateNormalFitResult(multivariateNormalMixtureExpectationMaximization.getFittedModel(), true, true);
}
 

import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
@Test
public void testRandomSelectFlatProbability() {
    final RandomGenerator rg = RandomGeneratorFactory.createRandomGenerator(new Random(13));
    final int NUM_SAMPLES = 1000;
    final List<Integer> choices = Arrays.asList(0,1,2);
    final List<Integer> result = IntStream.range(0, NUM_SAMPLES)
            .map(n -> GATKProtectedMathUtils.randomSelect(choices, j -> 1.0 / choices.size(), rg))
            .boxed()
            .collect(Collectors.toList());
    Assert.assertEquals(result.stream().filter(n -> n == 0).count(), NUM_SAMPLES / choices.size(), 50);
}
 

import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
@Test
public void testRandomSelect() {
    final RandomGenerator rg = RandomGeneratorFactory.createRandomGenerator(new Random(13));
    final int NUM_SAMPLES = 1000;
    final List<Integer> choices = Arrays.asList(-1,0,1);
    final List<Integer> result = IntStream.range(0, NUM_SAMPLES)
            .map(n -> GATKProtectedMathUtils.randomSelect(choices, j -> j*j/2.0, rg))
            .boxed()
            .collect(Collectors.toList());
    Assert.assertEquals(result.stream().filter(n -> n==0).count(), 0);
    Assert.assertEquals(result.stream().filter(n -> n == 1).count(), NUM_SAMPLES / 2, 50);
}
 

import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
@Test
public void testSegmentation() {
    final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(563));

    final List<Double> trueWeights = Arrays.asList(0.2, 0.5, 0.3);
    final List<Double> trueLog2CopyRatios = Arrays.asList(-2.0, 0.0, 1.4);
    final double trueMemoryLength = 1e5;
    final double trueStandardDeviation = 0.2;

    final CopyRatioHMM trueModel = new CopyRatioHMM(trueLog2CopyRatios, trueWeights,
            trueMemoryLength, trueStandardDeviation);

    final int chainLength = 10000;
    final List<SimpleInterval> positions = randomPositions("chr1", chainLength, rng, trueMemoryLength/4);
    final List<Integer> trueStates = trueModel.generateHiddenStateChain(positions);
    final List<Double> trueLog2CopyRatioSequence = trueStates.stream().map(n -> trueLog2CopyRatios.get(n)).collect(Collectors.toList());

    final List<Double> data = trueLog2CopyRatioSequence.stream()
            .map(cr -> generateData(trueStandardDeviation, cr, rng)).collect(Collectors.toList());

    final List<Target> targets = positions.stream().map(Target::new).collect(Collectors.toList());
    final ReadCountCollection rcc = new ReadCountCollection(targets, Arrays.asList("SAMPLE"), new Array2DRowRealMatrix(data.stream().mapToDouble(x->x).toArray()));
    final CopyRatioSegmenter segmenter = new CopyRatioSegmenter(10, rcc);
    final List<ModeledSegment> segments = segmenter.getModeledSegments();

    final double[] segmentCopyRatios = segments.stream()
            .flatMap(s -> Collections.nCopies((int) s.getTargetCount(), s.getSegmentMeanInLog2CRSpace()).stream())
            .mapToDouble(x -> x).toArray();

    final double averageCopyRatioError = IntStream.range(0, trueLog2CopyRatioSequence.size())
            .mapToDouble(n -> Math.abs(segmentCopyRatios[n] - trueLog2CopyRatioSequence.get(n)))
            .average().getAsDouble();

    Assert.assertEquals(averageCopyRatioError, 0, 0.025);
}
 

import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
@Test
public void testChromosomesOnDifferentSegments() {
    final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(563));
    final double[] trueLog2CopyRatios = new double[] {-2.0, 0.0, 1.7};
    final double trueMemoryLength = 1e5;

    final double trueStandardDeviation = 0.2;

    // randomly set positions
    final int chainLength = 100;
    final List<SimpleInterval> positions = randomPositions("chr1", chainLength, rng, trueMemoryLength/4);
    positions.addAll(randomPositions("chr2", chainLength, rng, trueMemoryLength/4));
    positions.addAll(randomPositions("chr3", chainLength, rng, trueMemoryLength/4));

    final int trueState = 2;    //fix everything to the same state 2

    final List<Double> data = new ArrayList<>();
    for (int n = 0; n < positions.size(); n++) {
        final double copyRatio = trueLog2CopyRatios[trueState];
        final double observed = generateData(trueStandardDeviation, copyRatio, rng);
        data.add(observed);
    }

    final List<Target> targets = positions.stream().map(Target::new).collect(Collectors.toList());
    final ReadCountCollection rcc = new ReadCountCollection(targets, Arrays.asList("SAMPLE"), new Array2DRowRealMatrix(data.stream().mapToDouble(x->x).toArray()));
    final CopyRatioSegmenter segmenter = new CopyRatioSegmenter(10, rcc);
    final List<ModeledSegment> segments = segmenter.getModeledSegments();

    //check that each chromosome has at least one segment
    final int numDifferentContigsInSegments = (int) segments.stream().map(ModeledSegment::getContig).distinct().count();
    Assert.assertEquals(numDifferentContigsInSegments, 3);
}
 

import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
@Test
public void testSegmentation() {
    final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(563));

    final List<Double> trueWeights = Arrays.asList(0.2, 0.5, 0.3);
    final List<Double> trueMinorAlleleFractions = Arrays.asList(0.12, 0.32, 0.5);
    final double trueMemoryLength = 1e5;

    final AlleleFractionGlobalParameters trueParams = new AlleleFractionGlobalParameters(1.0, 0.01, 0.01);

    final AlleleFractionHMM trueModel = new AlleleFractionHMM(trueMinorAlleleFractions, trueWeights,
            trueMemoryLength, AllelicPanelOfNormals.EMPTY_PON, trueParams);

    // randomly set positions
    final int chainLength = 10000;
    final List<SimpleInterval> positions = CopyRatioSegmenterUnitTest.randomPositions("chr1", chainLength, rng, trueMemoryLength/4);
    final List<Integer> trueStates = trueModel.generateHiddenStateChain(positions);
    final List<Double> truthMinorFractions = trueStates.stream().map(trueModel::getMinorAlleleFraction).collect(Collectors.toList());
    final AllelicCountCollection counts = generateCounts(truthMinorFractions, positions, rng, trueParams);

    final AlleleFractionSegmenter segmenter = new AlleleFractionSegmenter(10, counts, AllelicPanelOfNormals.EMPTY_PON);
    final List<ModeledSegment> segments = segmenter.getModeledSegments();
    final double[] segmentMinorFractions = segments.stream()
            .flatMap(s -> Collections.nCopies((int) s.getTargetCount(), s.getSegmentMean()).stream())
            .mapToDouble(x->x).toArray();

    final double averageMinorFractionError = IntStream.range(0, truthMinorFractions.size())
            .mapToDouble(n -> Math.abs(segmentMinorFractions[n] - truthMinorFractions.get(n)))
            .average().getAsDouble();

    Assert.assertEquals(averageMinorFractionError, 0, 0.01);
}
 

import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
private Object[][] getUnivariateGaussianTargetsWithDropout(final double sigma, final double dropoutRate) {
    Random rng = new Random(337);
    final RandomGenerator randomGenerator = RandomGeneratorFactory.createRandomGenerator(rng);
    NormalDistribution n = new NormalDistribution(randomGenerator, 1, sigma);
    final int numDataPoints = 10000;
    final int numEventPoints = 2000;

    // Randomly select dropoutRate of targets and reduce by 25%-75% (uniformly distributed)
    UniformRealDistribution uniformRealDistribution = new UniformRealDistribution(randomGenerator, 0, 1.0);
    final List<ReadCountRecord.SingleSampleRecord> targetList = new ArrayList<>();
    for (int i = 0; i < numDataPoints; i++){
        double coverage = n.sample() + (i < (numDataPoints - numEventPoints) ? 0.0 : 0.5);
        if (uniformRealDistribution.sample() < dropoutRate) {
            double multiplier = .25 + uniformRealDistribution.sample()/2;
            coverage = coverage * multiplier;
        }
        targetList.add(new ReadCountRecord.SingleSampleRecord(new Target("arbitrary_name", new SimpleInterval("chr1", 100 + 2*i, 101 + 2 * i)), coverage));
    }

    HashedListTargetCollection<ReadCountRecord.SingleSampleRecord> targets = new HashedListTargetCollection<>(targetList);

    List<ModeledSegment> segments = new ArrayList<>();
    segments.add(new ModeledSegment(new SimpleInterval("chr1", 100, 16050), 8000, 1));
    segments.add(new ModeledSegment(new SimpleInterval("chr1", 16100, 20200), 2000, 1.5));

    return new Object [] []{ {targets, segments}};
}
 

import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
default List<S> generateHiddenStateChain(final List<T> positions) {
    final RandomGenerator rg = RandomGeneratorFactory.createRandomGenerator(new Random(RANDOM_SEED_FOR_CHAIN_GENERATION));
    final List<S> hiddenStates = hiddenStates();
    final List<S> result = new ArrayList<>(positions.size());

    final S initialState = MathUtils.randomSelect(hiddenStates, s -> Math.exp(logPriorProbability(s, positions.get(0))), rg);
    result.add(initialState);

    IntStream.range(1, positions.size()).forEach(n ->
        result.add(MathUtils.randomSelect(hiddenStates,
                s -> Math.exp(logTransitionProbability(result.get(n-1), positions.get(n - 1), s, positions.get(n))), rg)));
    return result;
}
 

import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
@Test
public void testRandomSelectFlatProbability() {
    final RandomGenerator rg = RandomGeneratorFactory.createRandomGenerator(new Random(13));
    final int NUM_SAMPLES = 1000;
    final List<Integer> choices = Arrays.asList(0,1,2);
    final List<Integer> result = IntStream.range(0, NUM_SAMPLES)
            .map(n -> MathUtils.randomSelect(choices, j -> 1.0 / choices.size(), rg))
            .boxed()
            .collect(Collectors.toList());
    Assert.assertEquals(result.stream().filter(n -> n == 0).count(), NUM_SAMPLES / choices.size(), 50);
}
 

import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
@Test
public void testRandomSelect() {
    final RandomGenerator rg = RandomGeneratorFactory.createRandomGenerator(new Random(13));
    final int NUM_SAMPLES = 1000;
    final List<Integer> choices = Arrays.asList(-1,0,1);
    final List<Integer> result = IntStream.range(0, NUM_SAMPLES)
            .map(n -> MathUtils.randomSelect(choices, j -> j*j/2.0, rg))
            .boxed()
            .collect(Collectors.toList());
    Assert.assertEquals(result.stream().filter(n -> n==0).count(), 0);
    Assert.assertEquals(result.stream().filter(n -> n == 1).count(), NUM_SAMPLES / 2, 50);
}
 

import org.apache.commons.math3.random.RandomGeneratorFactory; //导入依赖的package包/类
@Test
public void testInitialization() {
    final double meanBias = 1.2;
    final double biasVariance = 0.04;
    final double outlierProbability = 0.02;
    final AlleleFractionGlobalParameters globalParameters = new AlleleFractionGlobalParameters(meanBias, biasVariance, outlierProbability);
    final int numSegments = 100;
    final double averageHetsPerSegment = 50.;
    final double averageDepth = 50.;
    final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(RANDOM_SEED));

    final SampleLocatableMetadata metadata = new SimpleSampleLocatableMetadata(
            "test-sample",
            new SAMSequenceDictionary(IntStream.range(0, numSegments)
                    .mapToObj(i -> new SAMSequenceRecord("chr" + i + 1, 1000))
                    .collect(Collectors.toList())));
    final AlleleFractionSimulatedData simulatedData = new AlleleFractionSimulatedData(
            metadata, globalParameters, numSegments, averageHetsPerSegment, averageDepth, rng);

    final AlleleFractionSegmentedData data = simulatedData.getData();
    final AlleleFractionState initializedState = new AlleleFractionInitializer(data).getInitializedState();

    Assert.assertEquals(initializedState.meanBias(), meanBias, ABSOLUTE_TOLERANCE);
    Assert.assertEquals(initializedState.biasVariance(), biasVariance, ABSOLUTE_TOLERANCE);
    Assert.assertEquals(initializedState.outlierProbability(), outlierProbability, ABSOLUTE_TOLERANCE);

    final double averageMinorFractionError = IntStream.range(0, numSegments)
            .mapToDouble(s -> Math.abs(initializedState.segmentMinorFraction(s) - simulatedData.getTrueState().segmentMinorFraction(s)))
            .average().getAsDouble();
    Assert.assertEquals(averageMinorFractionError, 0, ABSOLUTE_TOLERANCE);
}