Hi,
I have been trying to build a script in Beast that uses the Yang Codon
model with MCMC operators, but it's started to crash as I've added
more operators to the script and now it won't compile and run as
written. I know from the error message when it crashes says that the
bitflip operator for the localClock.changes is the problem, but I'm
not sure how to fix it and I don't know if I have all of the
appropriate operators for this model and method. I am posting the
executable portion of the script below; if you have any suggestions
please let me know. Thanks!
<!-- The unique patterns for codon positions
1 -->
<mergePatterns id="CP1.patterns">
<!-- The unique patterns from 1 to end every
3 -->
<!--
npatterns=7315
-->
<patterns from="1" every="3">
<alignment idref="alignment"/>
</patterns>
</mergePatterns>
<!-- The unique patterns for codon positions
2 -->
<mergePatterns id="CP2.patterns">
<!-- The unique patterns from 2 to end every
3 -->
<!--
npatterns=7332
-->
<patterns from="2" every="3">
<alignment idref="alignment"/>
</patterns>
</mergePatterns>
<!-- The unique patterns for codon positions
3 -->
<mergePatterns id="CP3.patterns">
<!-- The unique patterns from 3 to end every
3 -->
<!--
npatterns=7387
-->
<patterns from="3" every="3">
<alignment idref="alignment"/>
</patterns>
</mergePatterns>
<!-- A prior assumption that the population size has remained
constant -->
<!-- throughout the time spanned by the
genealogy. -->
<constantSize id="constant" units="substitutions">
<populationSize>
<parameter id="constant.popSize" value="0.7" lower="0.0"
upper="Infinity"/>
</populationSize>
</constantSize>
<!-- Generate a random starting tree under the coalescent
process -->
<coalescentTree id="startingTree" rootHeight="0.7">
<taxa idref="taxa"/>
<constantSize idref="constant"/>
</coalescentTree>
<!-- Generate a tree
model -->
<treeModel id="treeModel">
<coalescentTree idref="startingTree"/>
<rootHeight>
<parameter id="treeModel.rootHeight"/>
</rootHeight>
<nodeHeights internalNodes="true">
<parameter id="treeModel.internalNodeHeights"/>
</nodeHeights>
<nodeHeights internalNodes="true" rootNode="true">
<parameter id="treeModel.allInternalNodeHeights"/>
</nodeHeights>
</treeModel>
<!-- Generate a coalescent
likelihood -->
<coalescentLikelihood id="coalescent">
<model>
<constantSize idref="constant"/>
</model>
<populationTree>
<treeModel idref="treeModel"/>
</populationTree>
</coalescentLikelihood>
<!-- The random local clock model (Drummond & Suchard,
2010) -->
<randomLocalClockModel id="branchRates" ratesAreMultipliers="false">
<treeModel idref="treeModel"/>
<rates>
<parameter id="localClock.relativeRates"/>
</rates>
<rateIndicator>
<parameter id="localClock.changes"/>
</rateIndicator>
<clockRate>
<parameter id="clock.rate" value="1.0"/>
</clockRate>
</randomLocalClockModel>
<sumStatistic id="rateChanges" name="rateChangeCount"
elementwise="true">
<parameter idref="localClock.changes"/>
</sumStatistic>
<rateStatistic id="meanRate" name="meanRate" mode="mean"
internal="true" external="true">
<treeModel idref="treeModel"/>
<randomLocalClockModel idref="branchRates"/>
</rateStatistic>
<rateStatistic id="coefficientOfVariation"
name="coefficientOfVariation" mode="coefficientOfVariation"
internal="true" external="true">
<treeModel idref="treeModel"/>
<randomLocalClockModel idref="branchRates"/>
</rateStatistic>
<rateCovarianceStatistic id="covariance" name="covariance">
<treeModel idref="treeModel"/>
<randomLocalClockModel idref="branchRates"/>
</rateCovarianceStatistic>
<!-- The HKY substitution model (Hasegawa, Kishino & Yano,
1985) -->
<HKYModel id="CP1.hky">
<frequencies>
<frequencyModel dataType="nucleotide">
<frequencies>
<parameter id="CP1.frequencies" value="0.25 0.25 0.25 0.25"/>
</frequencies>
</frequencyModel>
</frequencies>
<kappa>
<parameter id="CP1.kappa" value="2.0" lower="0.0" upper="Infinity"/
</kappa>
</HKYModel>
<!-- The HKY substitution model (Hasegawa, Kishino & Yano,
1985) -->
<HKYModel id="CP2.hky">
<frequencies>
<frequencyModel dataType="nucleotide">
<frequencies>
<parameter id="CP2.frequencies" value="0.25 0.25 0.25 0.25"/>
</frequencies>
</frequencyModel>
</frequencies>
<kappa>
<parameter id="CP2.kappa" value="2.0" lower="0.0" upper="Infinity"/
</kappa>
</HKYModel>
<!-- The HKY substitution model (Hasegawa, Kishino & Yano,
1985) -->
<HKYModel id="CP3.hky">
<frequencies>
<frequencyModel dataType="nucleotide">
<frequencies>
<parameter id="CP3.frequencies" value="0.25 0.25 0.25 0.25"/>
</frequencies>
</frequencyModel>
</frequencies>
<kappa>
<parameter id="CP3.kappa" value="2.0" lower="0.0" upper="Infinity"/
</kappa>
</HKYModel>
<!-- site
model -->
<siteModel id="CP1.siteModel">
<substitutionModel>
<HKYModel idref="CP1.hky"/>
</substitutionModel>
<relativeRate>
<parameter id="CP1.mu" value="1.0" lower="0.0" upper="5.0"/>
</relativeRate>
<gammaShape gammaCategories="4">
<parameter id="CP1.alpha" value="0.5" lower="0.0" upper="1000.0"/>
</gammaShape>
</siteModel>
<!-- site
model -->
<siteModel id="CP2.siteModel">
<substitutionModel>
<HKYModel idref="CP2.hky"/>
</substitutionModel>
<relativeRate>
<parameter id="CP2.mu" value="1.0" lower="0.0" upper="5.0"/>
</relativeRate>
<gammaShape gammaCategories="4">
<parameter id="CP2.alpha" value="0.5" lower="0.0" upper="1000.0"/>
</gammaShape>
</siteModel>
<!-- site
model -->
<siteModel id="CP3.siteModel">
<substitutionModel>
<HKYModel idref="CP3.hky"/>
</substitutionModel>
<relativeRate>
<parameter id="CP3.mu" value="1.0" lower="0.0" upper="5.0"/>
</relativeRate>
<gammaShape gammaCategories="4">
<parameter id="CP3.alpha" value="0.5" lower="0.0" upper="1000.0"/>
</gammaShape>
</siteModel>
<compoundParameter id="allMus">
<parameter idref="CP1.mu"/>
<parameter idref="CP2.mu"/>
<parameter idref="CP3.mu"/>
</compoundParameter>
<treeLikelihood id="CP1.treeLikelihood" useAmbiguities="false">
<mergePatterns idref="CP1.patterns"/>
<treeModel idref="treeModel"/>
<siteModel idref="CP1.siteModel"/>
<randomLocalClockModel idref="branchRates"/>
</treeLikelihood>
<treeLikelihood id="CP2.treeLikelihood" useAmbiguities="false">
<mergePatterns idref="CP2.patterns"/>
<treeModel idref="treeModel"/>
<siteModel idref="CP2.siteModel"/>
<randomLocalClockModel idref="branchRates"/>
</treeLikelihood>
<treeLikelihood id="CP3.treeLikelihood" useAmbiguities="false">
<mergePatterns idref="CP3.patterns"/>
<treeModel idref="treeModel"/>
<siteModel idref="CP3.siteModel"/>
<randomLocalClockModel idref="branchRates"/>
</treeLikelihood>
<!-- Define
operators -->
<operators id="operators">
<scaleOperator scaleFactor="0.75" weight="0.1">
<parameter idref="CP1.kappa"/>
</scaleOperator>
<scaleOperator scaleFactor="0.75" weight="0.1">
<parameter idref="CP2.kappa"/>
</scaleOperator>
<scaleOperator scaleFactor="0.75" weight="0.1">
<parameter idref="CP3.kappa"/>
</scaleOperator>
<deltaExchange delta="0.01" weight="0.1">
<parameter idref="CP1.frequencies"/>
</deltaExchange>
<deltaExchange delta="0.01" weight="0.1">
<parameter idref="CP2.frequencies"/>
</deltaExchange>
<deltaExchange delta="0.01" weight="0.1">
<parameter idref="CP3.frequencies"/>
</deltaExchange>
<scaleOperator scaleFactor="0.75" weight="0.1">
<parameter idref="CP1.alpha"/>
</scaleOperator>
<scaleOperator scaleFactor="0.75" weight="0.1">
<parameter idref="CP2.alpha"/>
</scaleOperator>
<scaleOperator scaleFactor="0.75" weight="0.1">
<parameter idref="CP3.alpha"/>
</scaleOperator>
<deltaExchange delta="0.75" parameterWeights="8064 8064 8064"
weight="3">
<parameter idref="allMus"/>
</deltaExchange>
<scaleOperator scaleFactor="0.75" weight="15">
<parameter idref="localClock.relativeRates"/>
</scaleOperator>
<bitFlipOperator weight="15">
<parameter idref="localClock.changes"/>
</bitFlipOperator>
<subtreeSlide size="0.06999999999999999" gaussian="true"
weight="15">
<treeModel idref="treeModel"/>
</subtreeSlide>
<narrowExchange weight="15">
<treeModel idref="treeModel"/>
</narrowExchange>
<wideExchange weight="3">
<treeModel idref="treeModel"/>
</wideExchange>
<wilsonBalding weight="3">
<treeModel idref="treeModel"/>
</wilsonBalding>
<scaleOperator scaleFactor="0.75" weight="3">
<parameter idref="treeModel.rootHeight"/>
</scaleOperator>
<uniformOperator weight="30">
<parameter idref="treeModel.internalNodeHeights"/>
</uniformOperator>
<scaleOperator scaleFactor="0.75" weight="3">
<parameter idref="constant.popSize"/>
</scaleOperator>
<upDownOperator scaleFactor="0.75" weight="3">
<up>
</up>
<down>
<parameter idref="treeModel.allInternalNodeHeights"/>
</down>
</upDownOperator>
</operators>
<!-- Define
MCMC -->
<mcmc id="mcmc" chainLength="10000000" autoOptimize="true">
<posterior id="posterior">
<prior id="prior">
<logNormalPrior mean="1.0" stdev="1.25" offset="0.0"
meanInRealSpace="false">
<parameter idref="CP1.kappa"/>
</logNormalPrior>
<logNormalPrior mean="1.0" stdev="1.25" offset="0.0"
meanInRealSpace="false">
<parameter idref="CP2.kappa"/>
</logNormalPrior>
<logNormalPrior mean="1.0" stdev="1.25" offset="0.0"
meanInRealSpace="false">
<parameter idref="CP3.kappa"/>
</logNormalPrior>
<poissonPrior mean="0.6931471805599453" offset="0.0">
<statistic idref="rateChanges"/>
</poissonPrior>
<gammaPrior shape="0.5" scale="2.0" offset="0.0">
<parameter idref="localClock.relativeRates"/>
</gammaPrior>
<oneOnXPrior>
<parameter idref="constant.popSize"/>
</oneOnXPrior>
<coalescentLikelihood idref="coalescent"/>
</prior>
<likelihood id="likelihood">
<treeLikelihood idref="CP1.treeLikelihood"/>
<treeLikelihood idref="CP2.treeLikelihood"/>
<treeLikelihood idref="CP3.treeLikelihood"/>
</likelihood>
</posterior>
<operators idref="operators"/>
<!-- write log to
screen -->
<log id="screenLog" logEvery="1000">
<column label="Posterior" dp="4" width="12">
<posterior idref="posterior"/>
</column>
<column label="Prior" dp="4" width="12">
<prior idref="prior"/>
</column>
<column label="Likelihood" dp="4" width="12">
<likelihood idref="likelihood"/>
</column>
<column label="rootHeight" sf="6" width="12">
<parameter idref="treeModel.rootHeight"/>
</column>
<column label="clock.rate" sf="6" width="12">
<parameter idref="clock.rate"/>
</column>
</log>
<!-- write log to
file -->
<log id="fileLog" logEvery="1000" fileName="concat-dna.fasta.log"
overwrite="false">
<posterior idref="posterior"/>
<prior idref="prior"/>
<likelihood idref="likelihood"/>
<parameter idref="treeModel.rootHeight"/>
<parameter idref="constant.popSize"/>
<parameter idref="CP1.kappa"/>
<parameter idref="CP2.kappa"/>
<parameter idref="CP3.kappa"/>
<parameter idref="CP1.frequencies"/>
<parameter idref="CP2.frequencies"/>
<parameter idref="CP3.frequencies"/>
<parameter idref="CP1.alpha"/>
<parameter idref="CP2.alpha"/>
<parameter idref="CP3.alpha"/>
<compoundParameter idref="allMus"/>
<parameter idref="clock.rate"/>
<sumStatistic idref="rateChanges"/>
<rateStatistic idref="meanRate"/>
<rateStatistic idref="coefficientOfVariation"/>
<rateCovarianceStatistic idref="covariance"/>
<treeLikelihood idref="CP1.treeLikelihood"/>
<treeLikelihood idref="CP2.treeLikelihood"/>
<treeLikelihood idref="CP3.treeLikelihood"/>
<coalescentLikelihood idref="coalescent"/>
</log>
<!-- write tree log to
file -->
<logTree id="treeFileLog" logEvery="1000" nexusFormat="true"
fileName="concat-dna.fasta.trees" sortTranslationTable="true">
<treeModel idref="treeModel"/>
<randomLocalClockModel idref="branchRates"/>
<posterior idref="posterior"/>
</logTree>
</mcmc>
<report>
<property name="timer">
<mcmc idref="mcmc"/>
</property>
</report>
</beast>
