PoonLab / Kaphi

Kernel-embedded ABC-SMC for phylodynamic inference
GNU Affero General Public License v3.0
4 stars 2 forks source link

Drop particles where mu is greater than lambda (speciation models) #90

Closed MathiasRenaud closed 7 years ago

MathiasRenaud commented 7 years ago

run.smc stalls (but does not stop) when sample.priors returns a value of mu that is greater than the value of lambda for a particle.

First there will be a check that the given model is a speciation model. Then it will check if lambda and mu are both parameters before dropping particles where lambda is less than mu.

MathiasRenaud commented 7 years ago

Ideally, I'd like to set the weight of these invalid particles to 0, so that they would be ignored by .resample.particles and .perturb.particles. However, doing so causes an error in uniroot within next.epsilon before .perturb.particles can be reached.

This is because when the weight for a particle is 0, the value of epsilon.obj.func at both ends of the uniroot interval is negative.

I have a few other ideas for this issue which I am going to test out.

ArtPoon commented 7 years ago

We need a way to specify constraints on prior distributions in the yaml file. For example, we could do something like this:

priors:
  lambda:
    dist: 'gamma'
    hyperparameters:
    - shape: 2.0 # kappa
    - rate: 1.0  # theta
  mu:
    dist: 'gamma'
    hyperparameters:
    - shape: 2.0
    - rate: 1.0
   constraints:
     - "mu<lambda"  # this should get parsed

In order for this to work, we need to:

MathiasRenaud commented 7 years ago

The constraints heading is located in the yaml file between priors and proposals. The format is:

constraints:
  - "param1 < param2"

This string is parsed to the expression:

expression(theta['param1']<theta['param2'])

The expression is used in sample.priors and propose to recursively call the function if the constraint is not satisfied.

The parser splits the user's string at the whitespace so any order of operators and parameter names work.