use of dr.evomodel.tree.TreeHeightStatistic in project beast-mcmc by beast-dev.
the class RLYModelTest method randomLocalYuleTester.
private void randomLocalYuleTester(TreeModel treeModel, Parameter I, Parameter b, OperatorSchedule schedule) {
MCMC mcmc = new MCMC("mcmc1");
MCMCOptions options = new MCMCOptions(1000000);
TreeLengthStatistic tls = new TreeLengthStatistic(TL, treeModel);
TreeHeightStatistic rootHeight = new TreeHeightStatistic(TREE_HEIGHT, treeModel);
Parameter m = new Parameter.Default("m", 1.0, 0.0, Double.MAX_VALUE);
SpeciationModel speciationModel = new RandomLocalYuleModel(b, I, m, false, Units.Type.YEARS, 4);
Likelihood likelihood = new SpeciationLikelihood(treeModel, speciationModel, "randomYule.like");
ArrayLogFormatter formatter = new ArrayLogFormatter(false);
MCLogger[] loggers = new MCLogger[2];
loggers[0] = new MCLogger(formatter, 100, false);
loggers[0].add(likelihood);
loggers[0].add(rootHeight);
loggers[0].add(tls);
loggers[0].add(I);
loggers[1] = new MCLogger(new TabDelimitedFormatter(System.out), 100000, false);
loggers[1].add(likelihood);
loggers[1].add(rootHeight);
loggers[1].add(tls);
loggers[1].add(I);
mcmc.setShowOperatorAnalysis(true);
mcmc.init(options, likelihood, schedule, loggers);
mcmc.run();
List<Trace> traces = formatter.getTraces();
ArrayTraceList traceList = new ArrayTraceList("yuleModelTest", traces, 0);
for (int i = 1; i < traces.size(); i++) {
traceList.analyseTrace(i);
}
TraceCorrelation tlStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("root." + birthRateIndicator));
System.out.println("mean = " + tlStats.getMean());
System.out.println("expected mean = 0.5");
assertExpectation("root." + birthRateIndicator, tlStats, 0.5);
}
use of dr.evomodel.tree.TreeHeightStatistic in project beast-mcmc by beast-dev.
the class YuleModelTest method yuleTester.
// public void testYuleWithWideExchange() {
//
// TreeModel treeModel = new TreeModel("treeModel", tree);
// Doesn't compile...
// yuleTester(treeModel, ExchangeOperatorTest.getWideExchangeSchedule(treeModel));
// }
private void yuleTester(TreeModel treeModel, OperatorSchedule schedule) {
MCMC mcmc = new MCMC("mcmc1");
MCMCOptions options = new MCMCOptions(1000000);
TreeLengthStatistic tls = new TreeLengthStatistic(TL, treeModel);
TreeHeightStatistic rootHeight = new TreeHeightStatistic(TREE_HEIGHT, treeModel);
Parameter b = new Parameter.Default("b", 2.0, 0.0, Double.MAX_VALUE);
Parameter d = new Parameter.Default("d", 0.0, 0.0, Double.MAX_VALUE);
SpeciationModel speciationModel = new BirthDeathGernhard08Model(b, d, null, BirthDeathGernhard08Model.TreeType.TIMESONLY, Units.Type.YEARS);
Likelihood likelihood = new SpeciationLikelihood(treeModel, speciationModel, "yule.like");
ArrayLogFormatter formatter = new ArrayLogFormatter(false);
MCLogger[] loggers = new MCLogger[2];
loggers[0] = new MCLogger(formatter, 100, false);
loggers[0].add(likelihood);
loggers[0].add(rootHeight);
loggers[0].add(tls);
loggers[1] = new MCLogger(new TabDelimitedFormatter(System.out), 100000, false);
loggers[1].add(likelihood);
loggers[1].add(rootHeight);
loggers[1].add(tls);
mcmc.setShowOperatorAnalysis(true);
mcmc.init(options, likelihood, schedule, loggers);
mcmc.run();
List<Trace> traces = formatter.getTraces();
ArrayTraceList traceList = new ArrayTraceList("yuleModelTest", traces, 0);
for (int i = 1; i < traces.size(); i++) {
traceList.analyseTrace(i);
}
// expectation of root height for 4 tips and lambda = 2
// rootHeight = 0.541666
// TL = 1.5
TraceCorrelation tlStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TL));
assertExpectation(TL, tlStats, 1.5);
TraceCorrelation treeHeightStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TREE_HEIGHT));
assertExpectation(TREE_HEIGHT, treeHeightStats, 0.5416666);
}
use of dr.evomodel.tree.TreeHeightStatistic in project beast-mcmc by beast-dev.
the class OperatorAssert method irreducibilityTester.
private void irreducibilityTester(Tree tree, int numLabelledTopologies, int chainLength, int sampleTreeEvery) throws IOException, Importer.ImportException {
MCMC mcmc = new MCMC("mcmc1");
MCMCOptions options = new MCMCOptions(chainLength);
TreeModel treeModel = new TreeModel("treeModel", tree);
TreeLengthStatistic tls = new TreeLengthStatistic(TL, treeModel);
TreeHeightStatistic rootHeight = new TreeHeightStatistic(TREE_HEIGHT, treeModel);
OperatorSchedule schedule = getOperatorSchedule(treeModel);
Parameter b = new Parameter.Default("b", 2.0, 0.0, Double.MAX_VALUE);
Parameter d = new Parameter.Default("d", 0.0, 0.0, Double.MAX_VALUE);
SpeciationModel speciationModel = new BirthDeathGernhard08Model(b, d, null, BirthDeathGernhard08Model.TreeType.UNSCALED, Units.Type.YEARS);
Likelihood likelihood = new SpeciationLikelihood(treeModel, speciationModel, "yule.like");
MCLogger[] loggers = new MCLogger[2];
// loggers[0] = new MCLogger(new ArrayLogFormatter(false), 100, false);
// loggers[0].add(likelihood);
// loggers[0].add(rootHeight);
// loggers[0].add(tls);
loggers[0] = new MCLogger(new TabDelimitedFormatter(System.out), 10000, false);
loggers[0].add(likelihood);
loggers[0].add(rootHeight);
loggers[0].add(tls);
File file = new File("yule.trees");
file.deleteOnExit();
FileOutputStream out = new FileOutputStream(file);
loggers[1] = new TreeLogger(treeModel, new TabDelimitedFormatter(out), sampleTreeEvery, true, true, false);
mcmc.setShowOperatorAnalysis(true);
mcmc.init(options, likelihood, schedule, loggers);
mcmc.run();
out.flush();
out.close();
Set<String> uniqueTrees = new HashSet<String>();
HashMap<String, Integer> topologies = new HashMap<String, Integer>();
HashMap<String, HashMap<String, Integer>> treeCounts = new HashMap<String, HashMap<String, Integer>>();
NexusImporter importer = new NexusImporter(new FileReader(file));
int sampleSize = 0;
while (importer.hasTree()) {
sampleSize++;
Tree t = importer.importNextTree();
String uniqueNewick = TreeUtils.uniqueNewick(t, t.getRoot());
String topology = uniqueNewick.replaceAll("\\w+", "X");
if (!uniqueTrees.contains(uniqueNewick)) {
uniqueTrees.add(uniqueNewick);
}
HashMap<String, Integer> counts;
if (topologies.containsKey(topology)) {
topologies.put(topology, topologies.get(topology) + 1);
counts = treeCounts.get(topology);
} else {
topologies.put(topology, 1);
counts = new HashMap<String, Integer>();
treeCounts.put(topology, counts);
}
if (counts.containsKey(uniqueNewick)) {
counts.put(uniqueNewick, counts.get(uniqueNewick) + 1);
} else {
counts.put(uniqueNewick, 1);
}
}
TestCase.assertEquals(numLabelledTopologies, uniqueTrees.size());
TestCase.assertEquals(sampleSize, chainLength / sampleTreeEvery + 1);
Set<String> keys = topologies.keySet();
double ep = 1.0 / topologies.size();
for (String topology : keys) {
double ap = ((double) topologies.get(topology)) / (sampleSize);
// assertExpectation(ep, ap, sampleSize);
HashMap<String, Integer> counts = treeCounts.get(topology);
Set<String> trees = counts.keySet();
double MSE = 0;
double ep1 = 1.0 / counts.size();
for (String t : trees) {
double ap1 = ((double) counts.get(t)) / (topologies.get(topology));
// assertExpectation(ep1, ap1, topologies.get(topology));
MSE += (ep1 - ap1) * (ep1 - ap1);
}
MSE /= counts.size();
System.out.println("The Mean Square Error for the topolgy " + topology + " is " + MSE);
}
}
use of dr.evomodel.tree.TreeHeightStatistic in project beast-mcmc by beast-dev.
the class TreeHeightStatisticParser method parseXMLObject.
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
final String name = xo.getAttribute(Statistic.NAME, xo.getId());
final Tree tree = (Tree) xo.getChild(Tree.class);
return new TreeHeightStatistic(name, tree);
}
Aggregations