use of beast.evolution.tree.SCMigrationModel in project MultiTypeTree by tgvaughan.
the class TWB_TS_Test method test1.
@Test
public void test1() throws Exception {
System.out.println("TWB_test 1");
// Fix seed.
Randomizer.setSeed(42);
// Assemble initial MultiTypeTree
String newickStr = "((1[&deme=0]:1,2[&deme=0]:1)[&deme=0]:1," + "3[&deme=0]:2)[&deme=0]:0;";
MultiTypeTreeFromNewick mtTree = new MultiTypeTreeFromNewick();
mtTree.initByName("value", newickStr, "typeLabel", "deme");
// Assemble migration model:
RealParameter rateMatrix = new RealParameter("0.1 0.1");
RealParameter popSizes = new RealParameter("7.0 7.0");
SCMigrationModel migModel = new SCMigrationModel();
migModel.initByName("rateMatrix", rateMatrix, "popSizes", popSizes, "typeSet", new TypeSet("A", "B"));
// Assemble distribution:
StructuredCoalescentTreeDensity distribution = new StructuredCoalescentTreeDensity();
distribution.initByName("migrationModel", migModel, "multiTypeTree", mtTree);
// Set up state:
State state = new State();
state.initByName("stateNode", mtTree);
// Set up operator:
TypedWilsonBalding operatorTWB = new TypedWilsonBalding();
operatorTWB.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel, "alpha", 0.2);
Operator operatorMTTS = new MultiTypeTreeScale();
operatorMTTS.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel, "scaleFactor", 0.8, "useOldTreeScaler", false);
// Set up stat analysis logger:
MultiTypeTreeStatLogger logger = new MultiTypeTreeStatLogger();
logger.initByName("multiTypeTree", mtTree, "burninFrac", 0.2, "logEvery", 1000);
// Set up MCMC:
MCMC mcmc = new MCMC();
mcmc.initByName("chainLength", "1000000", "state", state, "distribution", distribution, "operator", operatorTWB, "operator", operatorMTTS, "logger", logger);
// Run MCMC:
mcmc.run();
System.out.format("height mean = %s\n", logger.getHeightMean());
System.out.format("height var = %s\n", logger.getHeightVar());
System.out.format("height ESS = %s\n", logger.getHeightESS());
// Direct simulation:
double[] heights = UtilMethods.getSimulatedHeights(migModel, new IntegerParameter("0 0 0"));
double simHeightMean = DiscreteStatistics.mean(heights);
double simHeightVar = DiscreteStatistics.variance(heights);
// Compare results with simulation results:
boolean withinTol = (logger.getHeightESS() > 400) && (Math.abs(logger.getHeightMean() - simHeightMean) < 1.0) && (Math.abs(logger.getHeightVar() - simHeightVar) < 30);
Assert.assertTrue(withinTol);
}
use of beast.evolution.tree.SCMigrationModel in project MultiTypeTree by tgvaughan.
the class TWB_TS_Test method testTWB2.
@Test
public void testTWB2() throws Exception {
System.out.println("TWB_test 2");
// Fix seed.
Randomizer.setSeed(42);
// Assemble initial MultiTypeTree
String newickStr = "((1[&deme=1]:1,2[&deme=0]:1)[&deme=0]:1," + "3[&deme=0]:2)[&deme=0]:0;";
MultiTypeTreeFromNewick mtTree = new MultiTypeTreeFromNewick();
mtTree.initByName("value", newickStr, "typeLabel", "deme");
// Assemble migration model:
RealParameter rateMatrix = new RealParameter("0.1 0.1");
RealParameter popSizes = new RealParameter("7.0 7.0");
SCMigrationModel migModel = new SCMigrationModel();
migModel.initByName("rateMatrix", rateMatrix, "popSizes", popSizes, "typeSet", new TypeSet("A", "B"));
// Assemble distribution:
StructuredCoalescentTreeDensity distribution = new StructuredCoalescentTreeDensity();
distribution.initByName("migrationModel", migModel, "multiTypeTree", mtTree);
// Set up state:
State state = new State();
state.initByName("stateNode", mtTree);
// Set up operator:
TypedWilsonBalding operatorTWB = new TypedWilsonBalding();
operatorTWB.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel, "alpha", 0.2);
Operator operatorMTTS = new MultiTypeTreeScale();
operatorMTTS.initByName("weight", 1.0, "multiTypeTree", mtTree, "migrationModel", migModel, "scaleFactor", 0.8, "useOldTreeScaler", false);
// Set up stat analysis logger:
MultiTypeTreeStatLogger logger = new MultiTypeTreeStatLogger();
logger.initByName("multiTypeTree", mtTree, "burninFrac", 0.1, "logEvery", 1000);
// Set up MCMC:
MCMC mcmc = new MCMC();
mcmc.initByName("chainLength", "1000000", "state", state, "distribution", distribution, "operator", operatorTWB, "logger", logger);
// Run MCMC:
mcmc.run();
System.out.format("height mean = %s\n", logger.getHeightMean());
System.out.format("height var = %s\n", logger.getHeightVar());
System.out.format("height ESS = %s\n", logger.getHeightESS());
// Direct simulation:
double[] heights = UtilMethods.getSimulatedHeights(migModel, new IntegerParameter("1 0 0"));
double simHeightMean = DiscreteStatistics.mean(heights);
double simHeightVar = DiscreteStatistics.variance(heights);
// Compare analysis results with truth:
boolean withinTol = (logger.getHeightESS() > 400) && (Math.abs(logger.getHeightMean() - simHeightMean) < 1.0) && (Math.abs(logger.getHeightVar() - simHeightVar) < 30);
Assert.assertTrue(withinTol);
}
use of beast.evolution.tree.SCMigrationModel in project MultiTypeTree by tgvaughan.
the class SCSimTest method test.
@Test
public void test() throws Exception {
System.out.println("SCSim test");
Randomizer.setSeed(42);
// Set up migration model.
RealParameter rateMatrix = new RealParameter();
rateMatrix.initByName("dimension", 2, "value", "0.1 0.1");
RealParameter popSizes = new RealParameter();
popSizes.initByName("value", "7.0 7.0");
SCMigrationModel migrationModel = new SCMigrationModel();
migrationModel.initByName("rateMatrix", rateMatrix, "popSizes", popSizes, "typeSet", new TypeSet("A", "B"));
// Specify leaf types:
IntegerParameter leafTypes = new IntegerParameter();
leafTypes.initByName("value", "0 0 0");
// Generate ensemble:
int reps = 100000;
double[] heights = new double[reps];
for (int i = 0; i < reps; i++) {
beast.evolution.tree.StructuredCoalescentMultiTypeTree sctree;
sctree = new beast.evolution.tree.StructuredCoalescentMultiTypeTree();
sctree.initByName("migrationModel", migrationModel, "leafTypes", leafTypes);
heights[i] = sctree.getRoot().getHeight();
}
double meanHeights = DiscreteStatistics.mean(heights);
double varHeights = DiscreteStatistics.variance(heights);
boolean withinTol = (Math.abs(meanHeights - 19.2) < 0.2) && (Math.abs(varHeights - 310) < 20);
Assert.assertTrue(withinTol);
}
Aggregations