use of dr.inference.model.Likelihood in project beast-mcmc by beast-dev.
the class ThreadedCompoundLikelihoodParser method parseXMLObject.
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
ThreadedCompoundLikelihood compoundLikelihood = new ThreadedCompoundLikelihood();
for (int i = 0; i < xo.getChildCount(); i++) {
if (xo.getChild(i) instanceof Likelihood) {
compoundLikelihood.addLikelihood((Likelihood) xo.getChild(i));
} else {
Object rogueElement = xo.getChild(i);
throw new XMLParseException("An element (" + rogueElement + ") which is not a likelihood has been added to a " + THREADED_COMPOUND_LIKELIHOOD + " element");
}
}
double weight = xo.getAttribute(WEIGHT, 0.0);
if (weight < 0)
throw new XMLParseException("Robust weight must be non-negative.");
compoundLikelihood.setWeightFactor(Math.exp(-weight));
return compoundLikelihood;
}
use of dr.inference.model.Likelihood in project beast-mcmc by beast-dev.
the class InvariantOperatorParser method parseXMLObject.
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
double weight = xo.getDoubleAttribute(MCMCOperator.WEIGHT);
Parameter parameter = (Parameter) xo.getChild(Parameter.class);
Likelihood likelihood = (Likelihood) xo.getChild(Likelihood.class);
return new InvariantOperator.Rotation(parameter, weight, likelihood);
}
use of dr.inference.model.Likelihood in project beast-mcmc by beast-dev.
the class DirtyLikelihoodOperatorParser method parseXMLObject.
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
double weight = xo.getDoubleAttribute(MCMCOperator.WEIGHT);
Likelihood likelihood = (Likelihood) xo.getChild(Likelihood.class);
return new DirtyLikelihoodOperator(likelihood, weight);
}
use of dr.inference.model.Likelihood 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);
}
use of dr.inference.model.Likelihood in project beast-mcmc by beast-dev.
the class TestCalibratedYuleModel method yuleTester.
private void yuleTester(TreeModel treeModel, OperatorSchedule schedule, Parameter brParameter, double S, int chainLength) throws IOException, TreeUtils.MissingTaxonException {
MCMC mcmc = new MCMC("mcmc1");
MCMCOptions options = new MCMCOptions(chainLength);
TreeLengthStatistic tls = new TreeLengthStatistic(TL, treeModel);
TreeHeightStatistic rootHeight = new TreeHeightStatistic(TREE_HEIGHT, treeModel);
SpeciationModel speciationModel = new BirthDeathGernhard08Model("yule", brParameter, null, null, BirthDeathGernhard08Model.TreeType.UNSCALED, Units.Type.SUBSTITUTIONS, false);
Likelihood speciationLikelihood = new SpeciationLikelihood(treeModel, speciationModel, "yule.like");
Taxa halfTaxa = new Taxa();
for (int i = 0; i < taxa.getTaxonCount() / 2; i++) {
halfTaxa.addTaxon(new Taxon("T" + Integer.toString(i)));
}
TMRCAStatistic tmrca = new TMRCAStatistic("tmrca(halfTaxa)", treeModel, halfTaxa, false, false);
DistributionLikelihood logNormalLikelihood = new DistributionLikelihood(new LogNormalDistribution(M, S), // meanInRealSpace="false"
0);
logNormalLikelihood.addData(tmrca);
MonophylyStatistic monophylyStatistic = new MonophylyStatistic("monophyly(halfTaxa)", treeModel, halfTaxa, null);
BooleanLikelihood booleanLikelihood = new BooleanLikelihood();
booleanLikelihood.addData(monophylyStatistic);
//CompoundLikelihood
List<Likelihood> likelihoods = new ArrayList<Likelihood>();
likelihoods.add(speciationLikelihood);
likelihoods.add(logNormalLikelihood);
likelihoods.add(booleanLikelihood);
Likelihood prior = new CompoundLikelihood(0, likelihoods);
prior.setId(CompoundLikelihoodParser.PRIOR);
ArrayLogFormatter logformatter = new ArrayLogFormatter(false);
MCLogger[] loggers = new MCLogger[1];
loggers[0] = new MCLogger(logformatter, (int) (options.getChainLength() / 10000), false);
loggers[0].add(speciationLikelihood);
loggers[0].add(rootHeight);
loggers[0].add(tmrca);
loggers[0].add(tls);
loggers[0].add(brParameter);
mcmc.setShowOperatorAnalysis(false);
mcmc.init(options, prior, schedule, loggers);
mcmc.run();
List<Trace> traces = logformatter.getTraces();
ArrayTraceList traceList = new ArrayTraceList("yuleModelTest", traces, 1000);
for (int i = 1; i < traces.size(); i++) {
traceList.analyseTrace(i);
}
NumberFormatter formatter = new NumberFormatter(8);
TraceCorrelation tlStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TL));
TraceCorrelation treeHeightStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("tmrca(halfTaxa)"));
// out.write("tmrcaHeight = \t");
out.write(formatter.format(treeHeightStats.getMean()));
out.write("\t");
double expectedNodeHeight = Math.pow(Math.E, (M + (Math.pow(S, 2) / 2)));
// out.write("expectation = \t");
out.write(formatter.format(expectedNodeHeight));
out.write("\t");
double error = Math.abs((treeHeightStats.getMean() - expectedNodeHeight) / expectedNodeHeight);
NumberFormat percentFormatter = NumberFormat.getNumberInstance();
percentFormatter.setMinimumFractionDigits(5);
percentFormatter.setMinimumFractionDigits(5);
// out.write("error = \t");
out.write(percentFormatter.format(error));
out.write("\t");
// out.write("tl.ess = \t");
out.write(Double.toString(tlStats.getESS()));
System.out.println("tmrcaHeight = " + formatter.format(treeHeightStats.getMean()) + "; expectation = " + formatter.format(expectedNodeHeight) + "; error = " + percentFormatter.format(error) + "; tl.ess = " + tlStats.getESS());
}
Aggregations