Skip to content

Commit

Permalink
When having sampled ancestors, the removal probability r doesn't affe…
Browse files Browse the repository at this point in the history 
…ct the sampling proportion s. (that is, s=psi/(psi+mu)).
  • Loading branch information
@denisekuehnert
denisekuehnert committed Apr 20, 2016
1 parent d94ff0d commit 8bf5e2d
Show file tree
Hide file tree
Showing 3 changed files with 19 additions and 25 deletions.
28 changes: 13 additions & 15 deletions src/beast/evolution/speciation/BirthDeathMigrationModel.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ public void initAndValidate() {

if (originBranch==null) throw new RuntimeException("Error: Origin specified but originBranch missing!");

updateOrigin(coltree.getRoot());
updateOrigin(coltree.getRoot());

if (!Boolean.valueOf(System.getProperty("beast.resume")) && treeInput.get().getRoot().getHeight() >= origin.get().getValue())
throw new RuntimeException("Error: origin("+T+") must be larger than tree height("+coltree.getRoot().getHeight()+")!");
Expand Down Expand Up @@ -197,8 +197,8 @@ void setupIntegrators(){ // set up ODE's and integrators
state = i/totalIntervals;

rho[i]= rhoChanges>0?
rhoSamplingChangeTimes.contains(times[i]) ? rhos[rhos.length > n ? (rhoChanges+1)*state+index(times[i%totalIntervals], rhoSamplingChangeTimes) : state] : 0.
: rhos[0];
rhoSamplingChangeTimes.contains(times[i]) ? rhos[rhos.length > n ? (rhoChanges+1)*state+index(times[i%totalIntervals], rhoSamplingChangeTimes) : state] : 0.
: rhos[0];
}
}

Expand Down Expand Up @@ -347,10 +347,10 @@ public double calculateTreeLogLikelihood(TreeInterface tree) {
collectTimes(T);
setRho();

if (updateRates() < 0 || (times[totalIntervals-1] > T)) {
if (updateRates() < 0 || (times[totalIntervals-1] > T)) {
logP = Double.NEGATIVE_INFINITY;
return logP;
}
}

double[] noSampleExistsProp = new double[n];

Expand Down Expand Up @@ -563,20 +563,18 @@ public void transformParameters(){
//
// death[i] = ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state] - psi[i];
//
if (!SAModel || removalAffectsSamplingProportion.get())
psi[i] = p[p.length > n ? (samplingChanges+1)*state+index(times[i%totalIntervals], samplingRateChangeTimes) : state]
* ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state] ;

if (!SAModel)
death[i] = ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state] - psi[i];
if (!SAModel) {
psi[i] = p[p.length > n ? (samplingChanges + 1) * state + index(times[i % totalIntervals], samplingRateChangeTimes) : state]
* ds[ds.length > n ? (deathChanges + 1) * state + index(times[i % totalIntervals], deathRateChangeTimes) : state];

death[i] = ds[ds.length > n ? (deathChanges + 1) * state + index(times[i % totalIntervals], deathRateChangeTimes) : state] - psi[i];
}
else {
r[i] = removalProbabilities[removalProbabilities.length > n ? (rChanges+1)*state+index(times[i%totalIntervals], rChangeTimes) : state];

if (!removalAffectsSamplingProportion.get())
psi[i] = p[p.length > n ? (samplingChanges+1)*state+index(times[i%totalIntervals], samplingRateChangeTimes) : state]
* ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state]
/ (1+(r[i]-1)*p[p.length > n ? (samplingChanges+1)*state+index(times[i%totalIntervals], samplingRateChangeTimes) : state]);
psi[i] = p[p.length > n ? (samplingChanges+1)*state+index(times[i%totalIntervals], samplingRateChangeTimes) : state]
* ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state]
/ (1+(r[i]-1)*p[p.length > n ? (samplingChanges+1)*state+index(times[i%totalIntervals], samplingRateChangeTimes) : state]);


death[i] = ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state] - psi[i]*r[i];
Expand Down
13 changes: 6 additions & 7 deletions src/beast/evolution/speciation/BirthDeathMigrationModelUncoloured.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -695,18 +695,17 @@ public void transformParameters(){
birth[i] = R[R.length > n ? (birthChanges+1)*state+index(times[i%totalIntervals], birthRateChangeTimes) : state]
* ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state] ;

if (!SAModel || removalAffectsSamplingProportion.get())
psi[i] = p[p.length > n ? (samplingChanges+1)*state+index(times[i%totalIntervals], samplingRateChangeTimes) : state]
* ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state] ;
if (!SAModel) {
psi[i] = p[p.length > n ? (samplingChanges + 1) * state + index(times[i % totalIntervals], samplingRateChangeTimes) : state]
* ds[ds.length > n ? (deathChanges + 1) * state + index(times[i % totalIntervals], deathRateChangeTimes) : state];

if (!SAModel)
death[i] = ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state] - psi[i];
death[i] = ds[ds.length > n ? (deathChanges + 1) * state + index(times[i % totalIntervals], deathRateChangeTimes) : state] - psi[i];
}

else {
r[i] = removalProbabilities[removalProbabilities.length > n ? (rChanges+1)*state+index(times[i%totalIntervals], rChangeTimes) : state];

if (!removalAffectsSamplingProportion.get())
psi[i] = p[p.length > n ? (samplingChanges+1)*state+index(times[i%totalIntervals], samplingRateChangeTimes) : state]
psi[i] = p[p.length > n ? (samplingChanges+1)*state+index(times[i%totalIntervals], samplingRateChangeTimes) : state]
* ds[ds.length > n ? (deathChanges+1)*state+index(times[i%totalIntervals], deathRateChangeTimes) : state]
/ (1+(r[i]-1)*p[p.length > n ? (samplingChanges+1)*state+index(times[i%totalIntervals], samplingRateChangeTimes) : state]);

Expand Down
3 changes: 0 additions & 3 deletions src/beast/evolution/speciation/PiecewiseBirthDeathSamplingDistribution.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,9 +119,6 @@ public abstract class PiecewiseBirthDeathSamplingDistribution extends SpeciesTre
public Input<Integer> stateNumber =
new Input<>("stateNumber", "The number of states or locations", Input.Validate.REQUIRED);

public Input<Boolean> removalAffectsSamplingProportion =
new Input<>("removalAffectsSamplingProportion", "In R0 param, is samplingProportion = samplingRate/(r*samplingRate+deathRate)? Default=true. (Alternative: samplingProportion = samplingRate/(samplingRate+deathRate)) ", true);

public Input<RealParameter> adjustTimesInput =
new Input<>("adjustTimes", "Origin of MASTER sims which has to be deducted from the change time arrays");
// <!-- HACK ALERT for reestimation from MASTER sims: adjustTimes is used to correct the forward changetimes such that they don't include orig-root (when we're not estimating the origin) -->
Expand Down

0 comments on commit 8bf5e2d

Please sign in to comment.