use of dr.evomodel.arg.coalescent.ARGUniformPrior in project beast-mcmc by beast-dev.
the class ARGAddRemoveOperatorTest method flatPriorTester.
private void flatPriorTester(ARGModel arg, int chainLength, int sampleTreeEvery, double nodeCountSetting, double rootHeightAlpha, double rootHeightBeta, int maxCount) throws IOException, Importer.ImportException {
MCMC mcmc = new MCMC("mcmc1");
MCMCOptions options = new MCMCOptions(chainLength);
// double nodeCountSetting = 2.0;
// double rootHeightAlpha = 100;
// double rootHeightBeta = 0.5;
OperatorSchedule schedule = getSchedule(arg);
ARGUniformPrior uniformPrior = new ARGUniformPrior(arg, maxCount, arg.getExternalNodeCount());
PoissonDistribution poisson = new PoissonDistribution(nodeCountSetting);
DistributionLikelihood nodeCountPrior = new DistributionLikelihood(poisson, 0.0);
ARGReassortmentNodeCountStatistic nodeCountStatistic = new ARGReassortmentNodeCountStatistic("nodeCount", arg);
nodeCountPrior.addData(nodeCountStatistic);
DistributionLikelihood rootPrior = new DistributionLikelihood(new GammaDistribution(rootHeightAlpha, rootHeightBeta), 0.0);
CompoundParameter rootHeight = (CompoundParameter) arg.createNodeHeightsParameter(true, false, false);
rootPrior.addData(rootHeight);
List<Likelihood> likelihoods = new ArrayList<Likelihood>();
likelihoods.add(uniformPrior);
likelihoods.add(rootPrior);
likelihoods.add(nodeCountPrior);
CompoundLikelihood compoundLikelihood = new CompoundLikelihood(1, likelihoods);
compoundLikelihood.setId("likelihood1");
MCLogger[] loggers = new MCLogger[3];
loggers[0] = new MCLogger(new TabDelimitedFormatter(System.out), 10000, false);
loggers[0].add(compoundLikelihood);
loggers[0].add(arg);
File file = new File("test.args");
file.deleteOnExit();
FileOutputStream out = new FileOutputStream(file);
loggers[1] = new ARGLogger(arg, new TabDelimitedFormatter(out), sampleTreeEvery, "test");
ArrayLogFormatter formatter = new ArrayLogFormatter(false);
loggers[2] = new MCLogger(formatter, sampleTreeEvery, false);
loggers[2].add(arg);
arg.getRootHeightParameter().setId("root");
loggers[2].add(arg.getRootHeightParameter());
mcmc.setShowOperatorAnalysis(true);
mcmc.init(options, compoundLikelihood, schedule, loggers);
mcmc.run();
out.flush();
out.close();
List<Trace> traces = formatter.getTraces();
// Set<String> uniqueTrees = new HashSet<String>();
//
// NexusImporter importer = new NexusImporter(new FileReader(file));
// while (importer.hasTree()) {
// Tree t = importer.importNextTree();
// uniqueTrees.add(Tree.Utils.uniqueNewick(t, t.getRoot()));
// }
//
// TestCase.assertEquals(numTopologies, uniqueTrees.size()); List<Trace> traces = formatter.getTraces();
ArrayTraceList traceList = new ArrayTraceList("ARGTest", traces, 0);
for (int i = 1; i < traces.size(); i++) {
traceList.analyseTrace(i);
}
TraceCorrelation nodeCountStats = traceList.getCorrelationStatistics(1);
TraceCorrelation rootHeightStats = traceList.getCorrelationStatistics(4);
assertExpectation("nodeCount", nodeCountStats, poisson.truncatedMean(maxCount));
assertExpectation(TreeModelParser.ROOT_HEIGHT, rootHeightStats, rootHeightAlpha * rootHeightBeta);
}
Aggregations