Examples with GammaSiteModel dr.oldevomodel.sitemodel.GammaSiteModel used on opensource projects

Example 6 with GammaSiteModel

use of dr.oldevomodel.sitemodel.GammaSiteModel in project beast-mcmc by beast-dev.

the class GeneralSubstitutionModelTest method testGeneralSubstitutionModel.

public void testGeneralSubstitutionModel() {
    // Sub model
    FrequencyModel freqModel = new FrequencyModel(dataType, alignment.getStateFrequencies());
    // dimension="5" value="1.0"
    Parameter ratesPara = new Parameter.Default(GeneralSubstitutionModelParser.RATES, 5, 1.0);
    // relativeTo="5"
    GeneralSubstitutionModel generalSubstitutionModel = new GeneralSubstitutionModel(dataType, freqModel, ratesPara, 4);
    //siteModel
    GammaSiteModel siteModel = new GammaSiteModel(generalSubstitutionModel);
    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, null, null, false, false, true, false, false);
    treeLikelihood.setId(TreeLikelihoodParser.TREE_LIKELIHOOD);
    // Operators
    OperatorSchedule schedule = new SimpleOperatorSchedule();
    MCMCOperator operator = new ScaleOperator(ratesPara, 0.5);
    operator.setWeight(1.0);
    schedule.addOperator(operator);
    Parameter rootHeight = treeModel.getRootHeightParameter();
    rootHeight.setId(TREE_HEIGHT);
    operator = new ScaleOperator(rootHeight, 0.5);
    operator.setWeight(1.0);
    schedule.addOperator(operator);
    Parameter internalHeights = treeModel.createNodeHeightsParameter(false, true, false);
    operator = new UniformOperator(internalHeights, 10.0);
    schedule.addOperator(operator);
    operator = new SubtreeSlideOperator(treeModel, 1, 1, true, false, false, false, CoercionMode.COERCION_ON);
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.NARROW, treeModel, 1.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.WIDE, treeModel, 1.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    operator = new WilsonBalding(treeModel, 1.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    // Log
    ArrayLogFormatter formatter = new ArrayLogFormatter(false);
    MCLogger[] loggers = new MCLogger[2];
    loggers[0] = new MCLogger(formatter, 1000, false);
    loggers[0].add(treeLikelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(ratesPara);
    loggers[1] = new MCLogger(new TabDelimitedFormatter(System.out), 100000, false);
    loggers[1].add(treeLikelihood);
    loggers[1].add(rootHeight);
    loggers[1].add(ratesPara);
    // MCMC
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(10000000);
    mcmc.setShowOperatorAnalysis(true);
    mcmc.init(options, treeLikelihood, schedule, loggers);
    mcmc.run();
    // time
    System.out.println(mcmc.getTimer().toString());
    // Tracer
    List<Trace> traces = formatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("GeneralSubstitutionModelTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    //      <expectation name="likelihood" value="-1815.75"/>
    //		<expectation name="treeModel.rootHeight" value="6.42048E-2"/>
    //		<expectation name="rateAC" value="6.08986E-2"/>
    TraceCorrelation likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TreeLikelihoodParser.TREE_LIKELIHOOD));
    assertExpectation(TreeLikelihoodParser.TREE_LIKELIHOOD, likelihoodStats, -1815.75);
    TraceCorrelation treeHeightStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TREE_HEIGHT));
    assertExpectation(TREE_HEIGHT, treeHeightStats, 0.0640787258170083);
    TraceCorrelation rateACStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(GeneralSubstitutionModelParser.RATES + "1"));
    assertExpectation(GeneralSubstitutionModelParser.RATES + "1", rateACStats, 0.061071756742081366);
}

Example 7 with GammaSiteModel

use of dr.oldevomodel.sitemodel.GammaSiteModel in project beast-mcmc by beast-dev.

the class LikelihoodTest method testLikelihoodGTR.

public void testLikelihoodGTR() {
    System.out.println("\nTest Likelihood using GTR:");
    // Sub model
    Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    Variable<Double> rateACValue = new Parameter.Default(GTRParser.A_TO_C, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
    Variable<Double> rateAGValue = new Parameter.Default(GTRParser.A_TO_G, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
    Variable<Double> rateATValue = new Parameter.Default(GTRParser.A_TO_T, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
    Variable<Double> rateCGValue = new Parameter.Default(GTRParser.C_TO_G, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
    Variable<Double> rateCTValue = new Parameter.Default(GTRParser.C_TO_T, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
    Variable<Double> rateGTValue = new Parameter.Default(GTRParser.G_TO_T, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
    GTR gtr = new GTR(rateACValue, rateAGValue, rateATValue, rateCGValue, rateCTValue, rateGTValue, f);
    //siteModel
    Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
    GammaSiteModel siteModel = new GammaSiteModel(gtr, mu, null, 4, null);
    //treeLikelihood
    SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);
    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, null, null, false, false, true, false, false);
    assertEquals("treeLikelihoodGTR", format.format(-1969.14584), format.format(treeLikelihood.getLogLikelihood()));
}

Example 8 with GammaSiteModel

use of dr.oldevomodel.sitemodel.GammaSiteModel in project beast-mcmc by beast-dev.

the class LikelihoodTest method testLikelihoodHKY85I.

public void testLikelihoodHKY85I() {
    System.out.println("\nTest Likelihood using HKY85I:");
    // Sub model
    Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 38.564672, 0, 100);
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);
    //siteModel
    Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
    Parameter invar = new Parameter.Default(GammaSiteModelParser.PROPORTION_INVARIANT, 0.701211, 0, 1.0);
    GammaSiteModel siteModel = new GammaSiteModel(hky, mu, null, 4, invar);
    //treeLikelihood
    SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);
    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, null, null, false, false, true, false, false);
    assertEquals("treeLikelihoodHKY85I", format.format(-1789.91240), format.format(treeLikelihood.getLogLikelihood()));
}

Example 9 with GammaSiteModel

use of dr.oldevomodel.sitemodel.GammaSiteModel in project beast-mcmc by beast-dev.

the class LikelihoodTest method testLikelihoodJC69.

public void testLikelihoodJC69() {
    System.out.println("\nTest Likelihood using JC69:");
    // Sub model
    Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 1.0, 0, 100);
    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, null, null, false, false, true, false, false);
    assertEquals("treeLikelihoodJC69", format.format(-1992.20564), format.format(treeLikelihood.getLogLikelihood()));
}

Example 10 with GammaSiteModel

use of dr.oldevomodel.sitemodel.GammaSiteModel in project beast-mcmc by beast-dev.

the class Coevolve method runModel2.

private void runModel2(PatternList patternList, PrintWriter pw, Tree tree, SubstitutionModel substModel, final Parameter betaParameter) {
    final Parameter muParameter = new Parameter.Default(1.0);
    muParameter.setId("mu");
    SiteModel siteModel = new GammaSiteModel(substModel, muParameter, null, 1, null);
    TreeModel treeModel = new TreeModel(tree);
    final TreeLikelihood treeLikelihood = new TreeLikelihood(patternList, treeModel, siteModel, null, null, false, false, true, false, false);
    treeLikelihood.setId("likelihood");
    MultivariateFunction function = new MultivariateFunction() {

        public double evaluate(double[] argument) {
            betaParameter.setParameterValue(0, argument[0]);
            betaParameter.setParameterValue(1, argument[1]);
            muParameter.setParameterValue(0, argument[2]);
            double lnL = -treeLikelihood.getLogLikelihood();
            // System.err.println("" + argument[0] + "\t" + argument[1] + "\t" + argument[2] + "\t" + lnL);
            return lnL;
        }

        public int getNumArguments() {
            return 3;
        }

        public double getLowerBound(int n) {
            return 0.0;
        }

        public double getUpperBound(int n) {
            return 100.0;
        }
    };
    MultivariateMinimum optimizer = new ConjugateGradientSearch();
    double lnL = optimizer.findMinimum(function, new double[] { 1.0, 1.0, 1.0 }, 6, 6);
    pw.write(betaParameter.getParameterValue(0) + "\t");
    pw.write(betaParameter.getParameterValue(1) + "\t");
    pw.write(muParameter.getParameterValue(0) + "\t");
    pw.write(lnL + "\n");
    pw.flush();
    System.out.println("" + betaParameter.getParameterValue(0) + "\t" + betaParameter.getParameterValue(1) + "\t" + muParameter.getParameterValue(0) + "\t" + lnL);
}