Search in sources :

Example 1 with FrequencyModel

use of dr.oldevomodel.substmodel.FrequencyModel in project beast-mcmc by beast-dev.

the class RecomboGen method main.

public static void main(String[] argv) {
    // Simulate sequences on this tree to generate sequences at the top of the
    // recombination process.
    Parameter kappa = new Parameter.Default(1, 2);
    Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
    FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, freqModel);
    // create site model
    SiteModel siteModel = new GammaSiteModel(hky);
    // create branch rate model
    Parameter rate = new Parameter.Default(1, 1.0E-4);
    BranchRateModel branchRateModel = new StrictClockBranchRates(rate);
//        RecomboGen recomboGen = new RecomboGen(8.0E-2, 100, 1000, 1.0E-4,
//                new String[] { "taxon1_0", "taxon2_0", "taxon3_10", "taxon4_10", "taxon5_20", "taxon6_20"},
//                new double[] { 0, 0, 10, 10, 20, 20} );
//
//        recomboGen.generate();
}
Also used : FrequencyModel(dr.oldevomodel.substmodel.FrequencyModel) GammaSiteModel(dr.oldevomodel.sitemodel.GammaSiteModel) HKY(dr.oldevomodel.substmodel.HKY) Parameter(dr.inference.model.Parameter) GammaSiteModel(dr.oldevomodel.sitemodel.GammaSiteModel) SiteModel(dr.oldevomodel.sitemodel.SiteModel)

Example 2 with FrequencyModel

use of dr.oldevomodel.substmodel.FrequencyModel in project beast-mcmc by beast-dev.

the class SequenceSimulator method getDefaultSiteModel.

/** generate simple site model, for testing purposes **/
static SiteModel getDefaultSiteModel() {
    Parameter kappa = new Parameter.Default(1, 2);
    Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);
    return new GammaSiteModel(hky);
}
Also used : FrequencyModel(dr.oldevomodel.substmodel.FrequencyModel) GammaSiteModel(dr.oldevomodel.sitemodel.GammaSiteModel) HKY(dr.oldevomodel.substmodel.HKY) Parameter(dr.inference.model.Parameter)

Example 3 with FrequencyModel

use of dr.oldevomodel.substmodel.FrequencyModel in project beast-mcmc by beast-dev.

the class StrictClockTest method testStrictClock.

public void testStrictClock() throws Exception {
    Parameter popSize = new Parameter.Default(ConstantPopulationModelParser.POPULATION_SIZE, 380.0, 0, 38000.0);
    ConstantPopulationModel constantModel = createRandomInitialTree(popSize);
    CoalescentLikelihood coalescent = new CoalescentLikelihood(treeModel, null, new ArrayList<TaxonList>(), constantModel);
    coalescent.setId("coalescent");
    // clock model
    Parameter rateParameter = new Parameter.Default(StrictClockBranchRates.RATE, 2.3E-5, 0, 100.0);
    StrictClockBranchRates branchRateModel = new StrictClockBranchRates(rateParameter);
    // Sub model
    Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 1.0, 0, 100.0);
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);
    //siteModel
    GammaSiteModel siteModel = new GammaSiteModel(hky);
    Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
    siteModel.setMutationRateParameter(mu);
    //treeLikelihood
    SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);
    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, branchRateModel, null, false, false, true, false, false);
    treeLikelihood.setId(TreeLikelihoodParser.TREE_LIKELIHOOD);
    // Operators
    OperatorSchedule schedule = new SimpleOperatorSchedule();
    MCMCOperator operator = new ScaleOperator(kappa, 0.75);
    operator.setWeight(1.0);
    schedule.addOperator(operator);
    operator = new ScaleOperator(rateParameter, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter allInternalHeights = treeModel.createNodeHeightsParameter(true, true, false);
    operator = new UpDownOperator(new Scalable[] { new Scalable.Default(rateParameter) }, new Scalable[] { new Scalable.Default(allInternalHeights) }, 0.75, 3.0, CoercionMode.COERCION_ON);
    schedule.addOperator(operator);
    operator = new ScaleOperator(popSize, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter rootHeight = treeModel.getRootHeightParameter();
    rootHeight.setId(TREE_HEIGHT);
    operator = new ScaleOperator(rootHeight, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter internalHeights = treeModel.createNodeHeightsParameter(false, true, false);
    operator = new UniformOperator(internalHeights, 30.0);
    schedule.addOperator(operator);
    operator = new SubtreeSlideOperator(treeModel, 15.0, 1.0, true, false, false, false, CoercionMode.COERCION_ON);
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.NARROW, treeModel, 15.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.WIDE, treeModel, 3.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    operator = new WilsonBalding(treeModel, 3.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    //CompoundLikelihood
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    likelihoods.add(coalescent);
    Likelihood prior = new CompoundLikelihood(0, likelihoods);
    prior.setId(CompoundLikelihoodParser.PRIOR);
    likelihoods.clear();
    likelihoods.add(treeLikelihood);
    Likelihood likelihood = new CompoundLikelihood(-1, likelihoods);
    likelihoods.clear();
    likelihoods.add(prior);
    likelihoods.add(likelihood);
    Likelihood posterior = new CompoundLikelihood(0, likelihoods);
    posterior.setId(CompoundLikelihoodParser.POSTERIOR);
    // Log
    ArrayLogFormatter formatter = new ArrayLogFormatter(false);
    MCLogger[] loggers = new MCLogger[2];
    loggers[0] = new MCLogger(formatter, 500, false);
    loggers[0].add(posterior);
    loggers[0].add(treeLikelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(rateParameter);
    loggers[0].add(popSize);
    loggers[0].add(kappa);
    loggers[0].add(coalescent);
    loggers[1] = new MCLogger(new TabDelimitedFormatter(System.out), 10000, false);
    loggers[1].add(posterior);
    loggers[1].add(treeLikelihood);
    loggers[1].add(rootHeight);
    loggers[1].add(rateParameter);
    loggers[1].add(coalescent);
    // MCMC
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(1000000);
    mcmc.setShowOperatorAnalysis(true);
    mcmc.init(options, posterior, schedule, loggers);
    mcmc.run();
    // time
    System.out.println(mcmc.getTimer().toString());
    // Tracer
    List<Trace> traces = formatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("RandomLocalClockTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    //        <expectation name="posterior" value="-3928.71"/>
    //        <expectation name="clock.rate" value="8.04835E-4"/>
    //        <expectation name="constant.popSize" value="37.3762"/>
    //        <expectation name="hky.kappa" value="18.2782"/>
    //        <expectation name="treeModel.rootHeight" value="69.0580"/>
    //        <expectation name="treeLikelihood" value="-3856.59"/>
    //        <expectation name="coalescent" value="-72.1285"/>
    TraceCorrelation likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(CompoundLikelihoodParser.POSTERIOR));
    assertExpectation(CompoundLikelihoodParser.POSTERIOR, likelihoodStats, -3928.71);
    likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TreeLikelihoodParser.TREE_LIKELIHOOD));
    assertExpectation(TreeLikelihoodParser.TREE_LIKELIHOOD, likelihoodStats, -3856.59);
    TraceCorrelation treeHeightStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TREE_HEIGHT));
    assertExpectation(TREE_HEIGHT, treeHeightStats, 69.0580);
    TraceCorrelation kappaStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(HKYParser.KAPPA));
    assertExpectation(HKYParser.KAPPA, kappaStats, 18.2782);
    TraceCorrelation rateStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(StrictClockBranchRates.RATE));
    assertExpectation(StrictClockBranchRates.RATE, rateStats, 8.04835E-4);
    TraceCorrelation popStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(ConstantPopulationModelParser.POPULATION_SIZE));
    assertExpectation(ConstantPopulationModelParser.POPULATION_SIZE, popStats, 37.3762);
    TraceCorrelation coalescentStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("coalescent"));
    assertExpectation("coalescent", coalescentStats, -72.1285);
}
Also used : FrequencyModel(dr.oldevomodel.substmodel.FrequencyModel) CompoundLikelihood(dr.inference.model.CompoundLikelihood) Likelihood(dr.inference.model.Likelihood) TreeLikelihood(dr.oldevomodel.treelikelihood.TreeLikelihood) CoalescentLikelihood(dr.evomodel.coalescent.CoalescentLikelihood) TreeLikelihood(dr.oldevomodel.treelikelihood.TreeLikelihood) ExchangeOperator(dr.evomodel.operators.ExchangeOperator) ArrayList(java.util.ArrayList) MCMC(dr.inference.mcmc.MCMC) SubtreeSlideOperator(dr.evomodel.operators.SubtreeSlideOperator) StrictClockBranchRates(dr.evomodel.branchratemodel.StrictClockBranchRates) CoalescentLikelihood(dr.evomodel.coalescent.CoalescentLikelihood) MCMCOptions(dr.inference.mcmc.MCMCOptions) ArrayLogFormatter(dr.inference.loggers.ArrayLogFormatter) WilsonBalding(dr.evomodel.operators.WilsonBalding) TraceCorrelation(dr.inference.trace.TraceCorrelation) SitePatterns(dr.evolution.alignment.SitePatterns) ConstantPopulationModel(dr.evomodel.coalescent.ConstantPopulationModel) TaxonList(dr.evolution.util.TaxonList) CompoundLikelihood(dr.inference.model.CompoundLikelihood) TabDelimitedFormatter(dr.inference.loggers.TabDelimitedFormatter) Trace(dr.inference.trace.Trace) GammaSiteModel(dr.oldevomodel.sitemodel.GammaSiteModel) ArrayTraceList(dr.inference.trace.ArrayTraceList) HKY(dr.oldevomodel.substmodel.HKY) Parameter(dr.inference.model.Parameter) MCLogger(dr.inference.loggers.MCLogger)

Example 4 with FrequencyModel

use of dr.oldevomodel.substmodel.FrequencyModel in project beast-mcmc by beast-dev.

the class BinaryCovarionModelTest method testTransitionProbabilitiesUnequalBaseFreqsEqualRates.

public void testTransitionProbabilitiesUnequalBaseFreqsEqualRates() {
    // with alpha == 1, the transition probability should be the same as binary jukes cantor
    alpha.setParameterValue(0, 1.0);
    switchingRate.setParameterValue(0, 1.0);
    frequencies.setParameterValue(0, 0.25);
    frequencies.setParameterValue(1, 0.75);
    FrequencyModel freqModel = new FrequencyModel(TwoStates.INSTANCE, frequencies);
    GeneralSubstitutionModel modelToCompare = new GeneralSubstitutionModel(TwoStates.INSTANCE, freqModel, null, 0);
    model.setupMatrix();
    double[] matrix = new double[16];
    double[] m = new double[4];
    double[] pi = model.getFrequencyModel().getFrequencies();
    for (double distance = 0.01; distance <= 1.005; distance += 0.01) {
        model.getTransitionProbabilities(distance, matrix);
        modelToCompare.getTransitionProbabilities(distance, m);
        double pChange = (matrix[1] + matrix[3]) * pi[0] + (matrix[4] + matrix[6]) * pi[1] + (matrix[9] + matrix[11]) * pi[2] + (matrix[12] + matrix[14]) * pi[3];
        double pChange2 = m[1] * frequencies.getParameterValue(0) + m[2] * frequencies.getParameterValue(1);
        System.out.println(distance + "\t" + pChange2 + "\t" + pChange);
        assertEquals(pChange2, pChange, 1e-14);
    }
}
Also used : FrequencyModel(dr.oldevomodel.substmodel.FrequencyModel) GeneralSubstitutionModel(dr.oldevomodel.substmodel.GeneralSubstitutionModel)

Example 5 with FrequencyModel

use of dr.oldevomodel.substmodel.FrequencyModel in project beast-mcmc by beast-dev.

the class BinaryCovarionModelTest method testTransitionProbabilitiesUnequalBaseFreqsUnequalRateFreqsEqualRates.

public void testTransitionProbabilitiesUnequalBaseFreqsUnequalRateFreqsEqualRates() {
    // with alpha == 1, the transition probability should be the same as binary jukes cantor
    alpha.setParameterValue(0, 1.0);
    switchingRate.setParameterValue(0, 1.0);
    frequencies.setParameterValue(0, 0.25);
    frequencies.setParameterValue(1, 0.75);
    hiddenFrequencies.setParameterValue(0, 0.1);
    hiddenFrequencies.setParameterValue(1, 0.9);
    FrequencyModel freqModel = new FrequencyModel(TwoStates.INSTANCE, frequencies);
    GeneralSubstitutionModel modelToCompare = new GeneralSubstitutionModel(TwoStates.INSTANCE, freqModel, null, 0);
    model.setupMatrix();
    double[] matrix = new double[16];
    double[] m = new double[4];
    double[] pi = model.getFrequencyModel().getFrequencies();
    for (double distance = 0.01; distance <= 1.005; distance += 0.01) {
        model.getTransitionProbabilities(distance, matrix);
        modelToCompare.getTransitionProbabilities(distance, m);
        double pChange = (matrix[1] + matrix[3]) * pi[0] + (matrix[4] + matrix[6]) * pi[1] + (matrix[9] + matrix[11]) * pi[2] + (matrix[12] + matrix[14]) * pi[3];
        double pChange2 = m[1] * frequencies.getParameterValue(0) + m[2] * frequencies.getParameterValue(1);
        //System.out.println(distance + "\t" + pChange2 + "\t" + pChange);
        assertEquals(pChange2, pChange, 1e-14);
    }
}
Also used : FrequencyModel(dr.oldevomodel.substmodel.FrequencyModel) GeneralSubstitutionModel(dr.oldevomodel.substmodel.GeneralSubstitutionModel)

Aggregations

FrequencyModel (dr.oldevomodel.substmodel.FrequencyModel)59 Parameter (dr.inference.model.Parameter)44 GammaSiteModel (dr.oldevomodel.sitemodel.GammaSiteModel)20 HKY (dr.oldevomodel.substmodel.HKY)18 SitePatterns (dr.evolution.alignment.SitePatterns)16 TreeLikelihood (dr.oldevomodel.treelikelihood.TreeLikelihood)16 DataType (dr.evolution.datatype.DataType)11 GeneralSubstitutionModel (dr.oldevomodel.substmodel.GeneralSubstitutionModel)6 ExchangeOperator (dr.evomodel.operators.ExchangeOperator)5 SubtreeSlideOperator (dr.evomodel.operators.SubtreeSlideOperator)5 WilsonBalding (dr.evomodel.operators.WilsonBalding)5 ArrayLogFormatter (dr.inference.loggers.ArrayLogFormatter)5 MCLogger (dr.inference.loggers.MCLogger)5 TabDelimitedFormatter (dr.inference.loggers.TabDelimitedFormatter)5 MCMC (dr.inference.mcmc.MCMC)5 MCMCOptions (dr.inference.mcmc.MCMCOptions)5 ArrayTraceList (dr.inference.trace.ArrayTraceList)5 Trace (dr.inference.trace.Trace)5 TraceCorrelation (dr.inference.trace.TraceCorrelation)5 GTR (dr.oldevomodel.substmodel.GTR)5