Skip to content

Commit

Permalink
debug BetaApproximationLikelihood
Browse files Browse the repository at this point in the history
  • Loading branch information
@rbouckaert
rbouckaert committed May 12, 2015
1 parent 8ab4d5d commit dbe5363
Show file tree
Hide file tree
Showing 2 changed files with 107 additions and 20 deletions.
93 changes: 93 additions & 0 deletions examples/testBetaApprox.xml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,93 @@
<!-- Simulation test 1: species = 4, lineages/species = 2, samples = 1000
u=v=1.0 alpha=2, beta=0.2 lambda=10
-->

<beast version='2.0' namespace='snap:beast.util:beast.core.util:beast.evolution'>

<map name='snapprior'>snap.likelihood.SnAPPrior</map>
<map name='snaptreelikelihood'>snap.likelihood.BetaApproximationLikelihood</map>

<!-- n = 1000 -->
<data spec='snap.Data' id='alignment' dataType='integerdata' statecount='3'>
<sequence taxon='A' totalcount='3'>
2,2,2,2,2,0,0,0,2,2,2,2,2,0,0,0,0,2,0,2,2,0,2,0,0,0,2,2,0,2,2,2,0,0,0,0,0,2,0,2,0,0,2,2,2,0,0,2,0,2,0,2,0,0,0,0,2,0,2,0,2,0,0,0,0,2,0,2,2,0,0,2,2,0,2,2,2,2,2,2,0,0,2,0,0,2,0,2,2,0,2,0,2,2,0,2,0,0,2,0,2,2,0,0,0,2,2,2,0,2,2,2,0,2,2,0,0,0,0,2,0,2,0,2,2,2,2,2,0,0,0,0,2,2,2,2,2,2,0,0,2,2,0,2,0,2,0,0,2,0,0,2,2,0,0,0,0,0,2,2,2,2,2,0,0,0,2,0,2,0,0,0,0,0,0,0,2,2,2,2,0,0,0,2,2,2,0,0,0,0,2,2,2,0,0,2,0,2,0,0,2,2,0,0,2,0,0,0,2,2,0,0,2,0,0,2,2,0,0,0,2,0,0,0,0,0,2,0,0,2,0,0,2,2,2,0,2,0,0,2,0,2,2,2,2,0,2,2,2,2,2,2,0,2,0,2,2,2,2,0,0,2,2,2,0,0,2,2,2,0,2,0,0,0,0,2,0,2,0,2,2,2,0,0,0,0,2,0,2,2,2,0,0,2,2,2,2,2,0,2,2,0,0,0,0,0,0,0,2,2,2,2,2,0,0,0,2,0,2,0,0,2,0,0,2,0,2,2,2,0,0,2,0,0,0,0,2,2,2,2,2,2,0,0,2,2,0,2,0,2,2,0,0,0,2,2,0,0,2,0,0,2,2,0,0,2,2,2,2,2,0,0,2,2,2,2,0,0,2,0,2,0,2,2,0,0,2,2,2,0,2,0,2,2,2,2,0,0,2,2,0,2,2,0,0,0,2,2,2,0,0,2,0,0,0,0,2,0,0,0,2,0,2,0,2,0,0,0,2,2,2,0,2,0,2,0,0,2,0,2,0,0,2,2,2,0,0,2,2,2,2,2,2,0,0,2,0,0,2,0,0,0,2,2,2,0,0,0,2,0,2,2,2,0,0,0,2,0,0,0,2,0,0,0,0,2,0,0,0,2,2,2,0,2,0,2,2,0,2,0,0,0,2,2,0,2,2,2,0,2,2,2,0,2,2,0,2,0,2,0,0,2,2,0,0,0,0,2,2,2,2,0,2,0,0,2,0,0,0,2,0,0,2,2,0,2,2,0,2,0,1,0,2,0,0,0,2,2,2,0,2,0,2,2,2,0,2,0,2,0,2,2,0,2,2,2,0,2,2,0,0,0,0,2,2,0,0,0,0,0,0,2,2,0,0,2,0,0,1,2,2,0,2,0,0,0,0,2,0,2,2,0,0,2,0,0,0,0,2,2,2,0,2,0,0,0,0,0,0,2,2,0,0,0,2,0,2,0,2,0,2,0,2,0,0,2,0,0,2,0,2,0,0,2,0,2,2,0,2,0,0,2,0,0,0,2,2,0,0,0,0,2,2,0,0,2,0,2,2,0,0,0,2,0,0,2,2,0,2,0,0,2,0,0,0,0,2,2,2,0,0,0,0,0,0,0,0,0,2,2,2,2,2,2,0,2,2,0,2,0,0,0,0,0,2,2,0,0,2,0,0,2,2,0,0,0,0,2,0,2,2,0,0,2,0,2,2,2,2,2,2,2,0,0,0,0,2,0,0,2,2,2,2,2,2,0,2,2,2,0,2,0,0,0,0,2,2,0,0,2,2,0,0,2,0,0,0,0,0,2,2,0,0,0,0,0,0,0,2,0,0,0,0,0,0,2,2,0,2,0,0,0,2,0,2,0,0,2,2,0,0,0,0,2,2,2,0,2,2,0,2,0,2,0,2,0,0,2,0,2,2,0,0,0,0,0,2,0,0,0,2,2,0,2,2,0,0,2,2,0,2,2,2,2,0,2,0,2,0,2,0,0,2,0,2,0,2,0,2,2,0,0,0,0,0,2,0,2,2,0,2,0,2,2,2,2,0,0,2,0,2,2,0,0,0,0,2,2,0,0,0,0,2,0,2,2,2,2,2,2,2,2,0,0,0,0,0,2,2,0,2,2,0,2,2,0,0,0,2,0,0,2,0,2,2,0,0,2,2,2,0,2,2,0,2,2,0,0,2,0,2,2,2,0,0,2,2,0,2,2,2,2,2,1,2,2,0,2,2,0,0,2,0,2,0,2,2,2,2,2,0,0,0,0,0,2,2,0,2,2,
</sequence>
<sequence taxon='B' totalcount='3'>
2,2,2,2,2,0,0,0,2,2,2,2,2,0,0,0,0,2,2,2,2,0,2,0,0,0,2,2,0,2,2,2,0,0,0,0,0,2,0,2,0,0,2,2,2,0,0,2,0,2,0,2,0,0,0,0,2,0,2,0,2,0,0,0,2,2,0,2,2,0,0,2,2,0,2,2,2,2,2,2,0,0,2,0,0,2,0,2,2,0,2,0,2,0,0,2,0,0,2,0,2,2,0,0,0,2,2,2,0,2,2,2,0,2,2,0,0,0,0,2,0,2,0,2,2,2,2,2,2,0,0,0,2,2,2,2,2,2,0,0,2,2,0,2,0,2,0,0,2,0,0,2,2,0,0,0,0,0,2,2,2,2,2,0,0,0,2,0,2,0,0,0,0,0,0,0,2,2,2,2,0,0,0,2,2,2,0,0,0,0,2,2,2,0,0,2,0,2,0,0,2,2,0,0,2,0,0,0,2,2,0,0,2,0,0,2,2,0,0,0,0,0,0,0,0,0,2,0,0,2,0,0,2,2,2,0,2,0,0,2,0,2,2,2,2,0,2,2,2,2,2,2,0,2,0,2,2,2,2,0,0,2,2,2,0,0,2,2,2,0,2,0,0,0,0,2,0,2,0,2,2,2,0,0,0,0,2,0,2,2,0,0,0,2,2,2,2,0,0,2,2,0,0,0,0,0,0,0,2,2,2,2,0,0,0,0,2,0,2,0,0,2,0,0,2,2,2,2,2,0,0,2,0,0,0,0,2,2,2,2,2,2,0,0,2,2,0,2,0,2,2,0,0,2,2,2,0,0,2,0,0,2,2,0,0,2,2,2,2,2,0,0,2,2,2,2,0,0,2,0,2,2,2,2,0,0,2,2,2,0,2,0,2,2,0,2,0,0,2,2,0,2,2,0,0,0,2,2,2,0,0,2,0,0,0,0,2,0,0,0,2,0,2,0,0,0,0,0,2,2,2,0,2,0,2,0,0,2,0,2,0,0,2,2,2,0,0,2,2,2,2,2,2,2,0,2,0,0,2,2,0,0,2,2,2,0,0,0,2,0,2,2,2,0,0,0,2,0,0,0,2,0,0,0,0,2,2,0,0,2,2,2,2,2,2,2,2,0,2,0,0,0,2,2,0,2,2,2,0,2,2,0,0,2,2,0,2,0,2,0,0,2,2,0,0,0,0,2,2,2,2,0,2,0,0,2,0,0,0,2,0,0,2,0,0,2,2,0,0,0,0,0,2,0,0,0,2,2,2,0,2,0,2,2,2,0,2,0,2,0,2,2,0,2,2,2,0,2,0,0,0,0,2,2,2,0,0,0,0,0,0,2,2,0,0,2,0,0,0,2,2,0,2,0,0,0,0,2,0,2,2,0,0,2,0,0,0,0,2,2,2,0,2,0,0,0,0,0,0,2,2,0,0,0,2,0,2,2,2,0,2,0,2,0,0,2,0,0,0,0,2,0,0,2,0,2,2,0,0,0,0,2,0,0,0,2,2,0,0,0,0,2,2,0,0,2,0,2,2,0,0,0,2,0,0,2,2,0,2,0,0,2,0,0,0,0,2,2,2,0,0,0,0,0,0,0,0,0,2,2,2,2,2,2,0,2,2,0,2,0,0,0,0,0,2,2,0,0,2,0,0,2,2,0,0,0,0,2,0,2,0,0,0,0,0,0,2,2,2,2,2,2,0,0,0,0,0,0,0,2,2,2,2,2,2,0,2,2,2,2,2,0,0,0,0,2,2,0,2,2,2,0,0,2,0,0,0,0,0,2,2,0,0,0,0,0,0,0,2,0,0,0,0,0,0,2,2,0,2,0,0,0,2,0,2,0,0,2,2,1,0,0,0,2,2,2,0,2,2,0,0,0,2,0,2,0,0,2,0,0,2,0,0,0,0,0,2,0,0,0,2,2,0,2,2,0,0,0,0,0,2,2,2,2,0,2,0,2,2,2,0,0,2,0,2,0,2,0,2,2,0,0,0,0,0,2,0,2,2,0,2,0,2,2,2,2,0,0,2,0,2,2,0,0,0,0,2,2,0,0,0,0,2,0,2,2,2,2,2,2,2,2,0,0,0,0,0,2,2,0,2,2,0,2,2,0,0,0,2,0,0,2,0,2,2,0,0,2,2,2,0,2,2,0,2,2,0,0,2,0,2,2,2,0,0,2,2,0,2,2,2,2,2,0,0,2,0,2,2,0,1,2,0,2,0,2,2,0,2,2,0,1,0,0,0,2,2,0,2,2,
</sequence>
<sequence taxon='C' totalcount='3'>
1,2,2,2,0,0,0,2,2,2,2,2,2,0,0,0,0,2,2,2,2,0,2,0,0,0,2,2,0,2,2,0,0,0,0,0,0,2,0,2,0,0,2,2,0,0,0,2,0,2,0,2,0,0,0,0,2,0,2,0,2,0,0,0,0,2,0,2,0,0,0,2,2,0,2,2,2,2,2,2,0,0,2,2,0,2,0,2,2,0,2,0,2,2,0,2,0,0,2,0,2,2,0,0,0,2,2,2,0,0,2,2,0,2,2,0,0,0,0,2,0,2,0,2,2,2,2,2,0,0,0,0,2,2,2,2,0,2,0,0,2,2,0,0,0,2,2,0,2,0,0,2,2,0,0,0,0,0,2,0,2,2,2,0,0,0,2,0,2,0,0,0,0,2,0,0,2,2,2,2,0,0,0,2,2,2,0,0,2,0,2,2,2,0,0,2,0,2,0,0,0,2,0,0,0,0,0,0,0,2,0,0,2,2,0,2,2,0,0,0,0,2,0,0,0,0,2,0,0,2,0,0,2,2,2,0,0,2,0,2,0,2,2,2,2,0,2,2,2,2,0,2,0,0,0,2,2,2,2,0,0,2,0,2,0,0,2,2,2,0,2,2,0,0,2,2,0,2,0,2,2,2,0,0,0,0,2,0,2,2,2,0,0,2,2,2,2,2,0,0,2,0,2,0,0,0,0,0,2,2,2,2,0,0,0,0,2,0,2,0,0,2,0,2,2,2,2,2,2,0,0,2,0,0,0,0,2,2,2,2,2,2,0,0,2,2,0,2,0,2,2,0,0,0,2,2,0,0,2,0,2,2,2,0,0,2,2,2,2,2,0,0,2,2,2,2,2,0,2,0,2,2,2,2,0,0,2,2,0,0,2,0,2,0,2,2,0,0,2,2,0,2,2,0,0,0,2,2,2,0,0,2,0,0,0,0,2,0,0,0,2,0,2,0,0,0,0,0,2,2,2,0,2,0,2,0,0,2,0,2,0,0,2,2,2,0,0,2,0,2,2,2,0,2,0,0,0,0,2,0,0,0,2,2,2,0,0,0,2,0,2,2,2,2,0,0,2,0,0,0,2,0,0,0,0,2,2,0,0,0,2,2,2,2,2,2,2,0,2,0,0,0,2,2,0,2,2,2,0,2,2,0,0,0,2,0,2,0,2,0,0,2,2,0,0,0,0,2,2,2,2,0,2,0,0,2,0,0,0,2,0,0,2,2,0,2,2,2,2,0,0,0,2,0,0,0,2,2,2,0,2,2,2,2,2,0,2,0,2,0,2,2,0,2,2,2,0,2,2,0,0,0,2,2,0,0,0,0,0,0,0,2,2,0,0,2,0,0,0,2,0,0,2,0,0,0,0,2,0,2,2,0,0,2,0,0,0,0,2,2,2,2,2,0,0,0,0,0,0,2,2,0,0,0,2,0,2,2,2,0,2,0,2,0,0,2,0,0,2,0,2,0,0,2,0,2,2,0,2,0,0,2,0,0,0,2,2,0,0,0,0,2,2,0,2,2,0,2,2,0,0,0,2,0,0,2,2,0,2,0,0,2,0,0,0,0,2,0,2,0,0,0,0,0,0,0,0,0,2,2,2,0,2,2,0,2,2,2,2,0,0,0,0,0,2,2,0,0,2,0,0,2,2,0,0,0,0,2,2,2,0,0,0,0,0,0,2,2,2,2,2,2,0,0,0,0,0,0,0,2,2,2,2,2,2,0,2,2,2,2,2,0,0,2,2,0,2,0,0,2,2,0,0,2,0,0,0,0,0,2,2,0,0,0,0,0,2,0,2,0,0,0,0,0,0,2,0,0,2,0,2,0,2,0,2,0,0,2,2,0,0,0,2,2,2,2,0,2,2,0,0,0,0,0,2,0,0,2,0,0,2,0,0,0,2,0,2,0,0,0,2,2,0,2,2,0,0,2,0,0,2,2,2,2,0,2,0,2,0,2,0,0,2,0,2,0,2,0,2,2,0,0,0,0,0,2,0,2,2,0,2,0,2,2,2,2,0,0,0,0,2,2,0,0,0,2,2,2,0,0,0,0,2,0,2,2,2,2,2,2,2,0,0,0,2,0,0,2,2,0,2,2,0,2,2,0,0,0,2,0,0,2,0,2,2,0,0,2,2,0,0,0,2,0,2,2,2,0,2,0,0,2,2,0,0,2,2,0,2,2,2,0,2,0,0,2,0,2,2,0,0,2,0,2,0,2,2,0,2,2,0,0,0,0,0,2,2,0,2,2,
</sequence>
<sequence taxon='D' totalcount='3'>
2,2,2,2,0,0,0,0,2,2,2,2,2,0,0,0,0,2,2,2,2,0,2,0,0,0,2,2,0,2,2,0,0,0,0,0,0,2,0,2,0,0,2,2,0,0,0,2,0,2,0,2,0,0,0,0,2,0,2,0,2,0,0,0,0,2,0,2,0,0,0,2,2,0,2,2,2,2,2,2,0,0,2,2,0,2,0,2,2,0,2,0,2,2,0,2,0,0,2,0,2,2,0,0,0,2,2,2,0,0,2,2,0,2,2,0,0,0,1,2,0,2,0,2,2,2,2,2,0,0,0,0,2,2,2,2,0,2,0,0,2,2,0,2,0,2,0,0,2,0,0,2,2,0,0,0,0,0,2,0,2,2,2,0,0,0,2,0,2,0,0,0,0,2,0,0,2,2,2,2,1,0,0,2,2,2,0,0,2,0,2,2,2,0,0,2,0,2,0,0,1,2,0,0,0,0,0,0,0,2,0,0,2,2,0,2,2,0,0,0,0,2,0,0,0,0,2,0,0,2,0,0,2,2,2,0,0,2,0,2,0,2,2,2,2,0,2,2,2,2,0,2,0,0,0,2,2,2,2,0,0,2,0,2,0,0,2,2,2,0,2,2,0,0,2,2,0,2,0,2,2,2,0,0,0,0,2,0,2,2,2,0,0,2,2,2,2,2,0,0,2,0,2,0,0,0,0,0,2,2,2,2,0,0,0,0,2,0,2,0,0,2,0,2,2,2,2,2,2,0,0,2,0,0,0,0,2,2,2,2,2,2,0,0,2,2,0,2,0,2,2,0,0,0,2,2,0,0,2,0,2,2,2,0,0,2,2,2,2,2,0,0,2,2,2,2,0,0,2,0,2,2,2,2,0,0,2,2,0,0,2,0,2,0,2,2,0,0,2,2,0,2,2,0,0,0,2,2,2,0,0,2,0,0,0,0,2,0,0,0,2,0,2,0,0,0,0,0,2,2,2,0,2,0,2,0,0,2,2,2,0,0,2,2,2,0,0,2,0,2,2,2,0,2,0,0,0,0,2,0,0,0,2,2,2,0,0,0,2,0,2,2,2,2,0,0,2,0,0,0,2,0,0,0,0,2,2,0,0,2,2,2,2,2,2,2,2,0,2,0,0,0,2,2,0,2,0,2,0,2,2,0,0,0,2,0,2,0,2,0,0,2,2,0,0,0,0,2,2,2,2,0,2,0,1,2,0,0,0,2,0,0,2,2,0,2,2,2,2,0,0,0,2,0,0,0,2,2,2,0,2,2,2,2,2,0,2,0,2,0,2,2,0,2,2,2,0,2,2,0,0,0,2,2,0,0,0,0,0,0,0,2,2,0,0,2,0,0,0,2,0,0,2,0,0,0,0,2,0,2,2,0,0,2,0,0,0,0,2,1,2,2,2,0,0,0,0,0,0,2,2,0,0,0,2,0,2,2,2,0,2,0,2,0,0,2,0,0,2,0,2,0,0,2,0,2,2,0,2,0,0,2,0,0,0,2,2,0,0,0,0,2,2,0,2,2,0,2,2,0,0,0,2,0,0,2,2,0,2,0,0,2,0,0,0,0,2,0,2,0,0,0,0,0,0,0,0,0,2,2,2,2,2,2,0,2,2,2,2,0,0,0,0,0,2,2,0,0,2,0,0,2,2,0,0,0,0,2,2,2,0,0,0,0,0,0,2,2,2,2,2,2,0,0,0,0,0,0,0,2,2,2,2,2,2,0,2,2,2,2,2,0,0,2,2,0,2,0,0,2,2,0,0,2,0,0,0,0,0,2,2,0,0,1,0,0,2,0,2,0,0,0,0,0,0,2,0,0,2,0,2,0,2,0,2,0,0,2,2,0,0,0,2,2,2,2,0,2,2,0,0,0,2,0,2,0,0,1,0,0,2,0,0,0,2,0,2,0,2,1,2,2,0,2,2,0,0,2,0,0,2,2,2,2,0,2,0,2,0,2,0,0,2,0,2,0,2,0,2,2,0,0,0,0,0,2,0,2,2,0,2,0,2,2,2,2,0,0,0,0,2,2,0,0,0,2,2,2,0,0,0,0,2,0,2,2,2,2,2,2,2,0,0,0,2,0,0,2,2,0,2,2,0,2,2,1,0,0,2,0,0,2,0,2,2,0,0,2,2,0,0,2,2,0,2,2,2,0,1,0,0,2,2,0,0,2,2,0,2,2,2,0,2,0,0,2,0,2,2,0,0,2,0,2,0,2,2,0,2,2,0,0,0,0,0,2,2,0,2,2,
</sequence>


</data>

<run id="mcmc" spec='beast.core.MCMC' chainLength="1000000" preBurnin="0">
<state storeEvery='1000'>
<parameter name='stateNode' id='coalescenceRate' dimension='7' value='10'/>

<!-- the following values may or may not be appropriate. Read the manual for further directions -->
<parameter name='stateNode' id='v' value='1.0' lower='0.0'/>
<parameter name='stateNode' id='u' value='1.0' lower='0.0'/>
<parameter name='stateNode' id='alpha' value='2' lower='0.0'/>
<parameter name='stateNode' id='beta' value='200' lower='0.0'/>
<parameter name='stateNode' id='lambda' value='10' lower='0.0'/>
<stateNode clusterType="upgma" estimate="true" id="tree" nodetype="snap.NodeData" spec="beast.util.ClusterTree" taxa='@alignment'/>
</state>

<distribution id="posterior" spec='CompoundDistribution'>
<snapprior name='distribution' id='prior' alpha='@alpha' beta='@beta' lambda='@lambda' tree='@tree' coalescenceRate='@coalescenceRate'/>
<snaptreelikelihood name='distribution' id='treeLikelihood' data='@alignment' tree='@tree' coalescenceRate='@coalescenceRate'/>
</distribution>

<!--
<stateDistribution idref='prior'/>
-->

<operator spec='operators.NodeSwapper' weight='0.5' tree='@tree'/>
<operator spec='operators.NodeBudger' weight='5' size='0.5' tree='@tree'/>
<operator spec='operators.ScaleOperator' scaleFactor='0.25' weight='1' tree='@tree'/>
<operator spec='operators.GammaMover' scale='0.5' weight='0.5' coalescenceRate='@coalescenceRate'/>
<operator spec='operators.RateMixer' scaleFactors='0.25' weight='1' tree='@tree' coalescenceRate='@coalescenceRate'/>

<!-- estimate u and v -->
<operator spec='operators.MutationMover' window='0.1' weight='0.25' u='@u' v='@v'/>


<logger logEvery="1000">
<log idref="u"/>
<log idref="v"/>
<log idref="prior"/>
<log idref="treeLikelihood"/>
<log idref="posterior"/>
<log spec='snap.ThetaLogger' coalescenceRate='@coalescenceRate'/>
<log spec='beast.evolution.tree.TreeHeightLogger' tree='@tree'/>
</logger>
<logger logEvery="1000" fileName="testBetaApproximation.log">
<model idref='posterior'/>
<log idref="u"/>
<log idref="v"/>
<log idref="prior"/>
<log idref="treeLikelihood"/>
<log idref="posterior"/>
<log idref='coalescenceRate'/>
<log spec='snap.ThetaLogger' coalescenceRate='@coalescenceRate'/>
<log spec='beast.evolution.tree.TreeHeightLogger' tree='@tree'/>
<log spec='TreeLengthLogger' tree='@tree'/>
</logger>
<logger logEvery="1000" fileName="testBetaApproximation.trees" mode='tree'>
<log spec='beast.evolution.tree.TreeWithMetaDataLogger' tree="@tree" metadata='@coalescenceRate'/>
</logger>

<!--
<logger logEvery="1000" fileName="test2.$(seed).trees" mode='tree'>
<log spec='snap.CoalescentUnitTreeLogger' tree="@tree" coalescenceRate='@coalescenceRate'/>
</logger>
-->
</run>

</beast>

34 changes: 14 additions & 20 deletions src/snap/likelihood/BetaApproximationLikelihood.java
View file
Original file line number Diff line number Diff line change
@@ -1,32 +1,31 @@
package snap.likelihood;



import snap.Data;
import beast.core.Citation;
import beast.core.Description;
import beast.evolution.alignment.Alignment;
import beast.evolution.branchratemodel.BranchRateModel;
import beast.evolution.branchratemodel.StrictClockModel;
import beast.core.Input;
import beast.core.Input.Validate;
import beast.core.parameter.RealParameter;
import beast.evolution.likelihood.GenericTreeLikelihood;
import beast.evolution.sitemodel.SiteModel;
import beast.evolution.tree.Node;
import beast.evolution.tree.Tree;
import beast.math.Binomial;
import beast.math.GammaFunction;

import org.apache.commons.math.distribution.BetaDistribution;

@Description("Implementation of Siren et al. Beta SNP approximation using Hiscott et al. integrator")
@Citation("")
public class BetaApproximationLikelihood extends GenericTreeLikelihood {
public Input<RealParameter> coalescenceRateInput = new Input<RealParameter>("coalescenceRate", "population size parameter with one value for each node in the tree");

public BetaApproximationLikelihood() {
siteModelInput.setRule(Validate.OPTIONAL);
}

Tree tree;
<<<<<<< HEAD
Double rootAlpha; //Parameter for root distribution.
=======
Double rootTheta; //Parameter for root distribution.
>>>>>>> 998524460c43138ca5124e4ba57b30d99803bf95
Data data;

int nPoints; //Size of mesh used in integration
Expand Down Expand Up @@ -55,21 +54,20 @@ public class BetaApproximationLikelihood extends GenericTreeLikelihood {
*/
protected double[] m_branchLengths;
protected double[] storedBranchLengths;

RealParameter coalescenceRate;


@Override
public void initAndValidate() throws Exception {
tree = (Tree) treeInput.get();
data = dataInput.get();
data = (Data) dataInput.get();
patternLogLikelihoods = new double[data.getPatternCount()];
coalescenceRate = coalescenceRateInput.get();

nPoints = 20; //TODO: read this from the xml
nPoints = data.getTaxonCount();//20; //TODO: read this from the xml
int nNodes = tree.getNodeCount();
partialIntegral = new double[nPoints][nNodes];
<<<<<<< HEAD
=======
rootTheta = 0.5; //We should get this from the PARAMETERS.
>>>>>>> 998524460c43138ca5124e4ba57b30d99803bf95

hasDirt = Tree.IS_FILTHY;
}
Expand All @@ -82,7 +80,7 @@ public void initAndValidate() throws Exception {
*/
@Override
public double calculateLogP() {
rootTheta = coalescenceRate.getValue(tree.getRoot().getNr()); //We should get this from the PARAMETERS.

double logL=0.0;

Expand Down Expand Up @@ -239,11 +237,7 @@ private double transProb(double x,double y, double t) {
}

private double rootProb(double x) {
<<<<<<< HEAD
return betaDensity(x,rootAlpha,rootAlpha);
=======
return betaDensity(x,rootTheta,rootTheta);
>>>>>>> 998524460c43138ca5124e4ba57b30d99803bf95
}


Expand Down

0 comments on commit dbe5363

Please sign in to comment.