use of beast.core.parameter.RealParameter in project bacter by tgvaughan.
the class ACGLikelihoodTest method testLikelihoodUsingSimulatedData.
@Test
public void testLikelihoodUsingSimulatedData() throws Exception {
ConstantPopulation popFunc = new ConstantPopulation();
popFunc.initByName("popSize", new RealParameter("1.0"));
Locus locus = new Locus("locus", 10000);
TaxonSet taxonSet = getTaxonSet(10);
ConversionGraph acg = new SimulatedACG();
acg.initByName("rho", 5.0 / locus.getSiteCount(), "delta", 1000.0, "populationModel", popFunc, "locus", locus, "taxonset", taxonSet);
System.out.println(acg);
// Site model:
JukesCantor jc = new JukesCantor();
jc.initByName();
SiteModel siteModel = new SiteModel();
siteModel.initByName("mutationRate", new RealParameter("1"), "substModel", jc);
// Simulate alignment:
SimulatedAlignment alignment = new SimulatedAlignment();
alignment.initByName("acg", acg, "siteModel", siteModel, "outputFileName", "simulated_alignment.nexus", "useNexus", true);
// Calculate likelihood:
ACGLikelihood argLikelihood = new ACGLikelihood();
argLikelihood.initByName("locus", locus, "data", alignment, "tree", acg, "siteModel", siteModel);
double logP = argLikelihood.calculateLogP();
// Compare product of likelihoods of "marginal alignments" with
// likelihood computed using RGL.
ACGLikelihoodSlow argLikelihoodSlow = new ACGLikelihoodSlow();
argLikelihoodSlow.initByName("locus", locus, "data", alignment, "tree", acg, "siteModel", siteModel);
double logPprime = argLikelihoodSlow.calculateLogP();
double relError = 2.0 * Math.abs(logP - logPprime) / Math.abs(logP + logPprime);
System.out.format("logP=%g\nlogPprime=%g\nrelError=%g\n", logP, logPprime, relError);
assertTrue(relError < 1e-13);
}
use of beast.core.parameter.RealParameter in project bacter by tgvaughan.
the class SkylinePopulationFunctionTest method test5.
@Test
public void test5() throws Exception {
SkylinePopulationFunction skyline = new SkylinePopulationFunction();
skyline.initByName("acg", acg, "popSizes", new RealParameter("1.0 1.0 5.0 1.0"), "groupSizes", new IntegerParameter("0"), "nGridPoints", 5, "piecewiseLinear", true);
for (double t = 0.0; t < 10; t += 0.01) assertTrue(Math.abs(t - skyline.getInverseIntensity(skyline.getIntensity(t))) < 1e-14);
}
use of beast.core.parameter.RealParameter in project bacter by tgvaughan.
the class SkylinePopulationFunctionTest method test6.
@Test
public void test6() throws Exception {
SkylinePopulationFunction skyline = new SkylinePopulationFunction();
skyline.initByName("acg", acg, "popSizes", new RealParameter("1.0 1.0 5.0 1.0"), "groupSizes", new IntegerParameter("0"), "nGridPoints", 5, "piecewiseLinear", true);
for (CFEventList.Event cfEvent : acg.getCFEvents()) {
double t = cfEvent.getHeight();
assertTrue(Math.abs(t - skyline.getInverseIntensity(skyline.getIntensity(t))) < 1e-14);
}
}
use of beast.core.parameter.RealParameter in project bacter by tgvaughan.
the class AddRemoveConversionTest method testProbability.
/**
* Tests whether probability of proposing a conversion lines up with
* conversion probability found in ACGCoalescent.
* @throws java.lang.Exception
*/
@Test
public void testProbability() throws Exception {
ConstantPopulation popFunc = new ConstantPopulation();
popFunc.initByName("popSize", new RealParameter("1.0"));
Locus locus = new Locus("locus", 10000);
TaxonSet taxonSet = getTaxonSet(10);
SimulatedACG acg = new SimulatedACG();
acg.initByName("rho", 1.0 / locus.getSiteCount(), "delta", 50.0, "locus", locus, "taxonset", taxonSet, "populationModel", popFunc);
RealParameter rho = new RealParameter(Double.toString(1.0 / locus.getSiteCount()));
RealParameter delta = new RealParameter("50.0");
AddRemoveConversion operator = new AddRemoveConversion();
operator.initByName("weight", 1.0, "acg", acg, "delta", delta, "populationModel", popFunc);
ACGCoalescent coal = new ACGCoalescent();
coal.initByName("tree", acg, "populationModel", popFunc, "rho", rho, "delta", delta);
double logP1 = 0.0;
double logP2 = 0.0;
for (Conversion conv : acg.getConversions(locus)) {
logP1 += operator.getConversionProb(conv);
logP2 += coal.calculateConversionLogP(conv);
}
System.out.println("logP1 = " + logP1);
System.out.println("logP2 = " + logP2);
assertTrue(Math.abs(logP1 - logP2) < 1e-10);
}
use of beast.core.parameter.RealParameter in project bacter by tgvaughan.
the class PiecewisePopulationFunction method main.
/**
* Main method for testing.
*
* @param args command line arguments (unused)
*/
public static void main(String[] args) throws Exception {
String acgString = "[&15,0,1.3905355989030808,31,770,1.597708055397074] " + "[&30,931,2.4351280458424904,36,2486,3.78055549386568] " + "[&15,941,2.0439957300083322,38,2364,6.911056700367016] " + "[&36,1091,4.285505683622974,38,2589,9.867725913197855] " + "((((10:0.5385300170206817,(17:0.116794353049212," + "((3:0.039229346597297564,12:0.039229346597297564)23:0.04582913870888949," + "13:0.08505848530618705)24:0.03173586774302495)26:0.4217356639714697)28:1.8114199763246093," + "((8:0.10883006062265468,2:0.10883006062265468)25:0.556428062025291," + "(6:0.5393311342677402,11:0.5393311342677402)29:0.12592698838020555)31:1.6846918706973453)34:1.4536824928125807," + "(1:0.47184545557390367,14:0.47184545557390367)27:3.331787030583968)37:2.9704369411362554," + "(((15:2.0624287390593707,((16:0.01825347077733299,19:0.01825347077733299)21:0.7668749128372041," + "(7:0.008018731329538273,9:0.008018731329538273)20:0.7771096522849988)32:1.2773003554448337)33:0.7487092404613747," + "4:2.8111379795207454)35:0.1331794525400949,((0:0.0243537216663141," + "5:0.0243537216663141)22:0.5681537100482162,18:0.5925074317145304)30:2.35181000034631)36:3.829751995233287)38:0.0";
ConversionGraph acg = new ConversionGraph();
acg.initByName("siteCount", 10000, "fromString", acgString);
PiecewisePopulationFunction skyline = new PiecewisePopulationFunction();
skyline.initByName("acg", acg, "popSizes", new RealParameter("1.0 1.0 5.0 1.0 2.0"), "groupSizes", new IntegerParameter("0"), "piecewiseLinear", true);
try (PrintStream ps = new PrintStream("out.txt")) {
ps.println("t N intensity intensityInv");
double t = 0.0;
while (t < 10) {
ps.format("%g %g %g %g\n", t, skyline.getPopSize(t), skyline.getIntensity(t), skyline.getInverseIntensity(skyline.getIntensity(t)));
t += 0.001;
}
ps.close();
}
}
Aggregations