Distribution for waiting until an event?

[bredelings]

The new front page for https://treeppl.org/ looks very nice. The type system and the ability to walk trees recursively seems very attractive compared to (say) revbayes.

The separation of function from model function also seems nice. It seems to mirror the fact that (in Haskell terminology) these two functions operate in different monads. I wonder if adding a comment about this to the code on the from page would make sense, because at first I did not see the model keyword.

Hmm... on second thought, the fact that the "observe" function can run in non-model function functions seems slightly odd. I presume that this results from the fact that no particles are created by observing, only by sampling? For MCMC, I guess sampling creates new random variables in the MCMC state, whereas observing does not (or does not need to). I don't know.

The main weird-looking thing about this code is observe 0 and observe 0.0 . It took some thinking to understand what this is supposed to mean. And even now I'm not totally sure this is correct: are all the tips at the present -- I mean, is node.age == 0.0 for all tips? If not, then where do these tip ages come from? It looks like there is no extinction, and no random observation of ancestors. So, the node times would have to be fixed in advance. But then why do they appear to cause extinction? I'm not sure I follow this.

Does TreePPL have an optional type similar to Haskell's Maybe a or Rusts's Option<a>? If so, perhaps you could have a distribution PoissonFirstEvent(lambda, duration) that returns either Some(time) or None depending on whether the event happens before the end of the interval or not. Then you could have something like

function walk(node: Tree, age: Real, lambda: Real) {
  # in a constant-rate birth process there are no extinctions, so all tips are at time 0.
  if node is Node {
    # a speciation happened before the present
    observe Some(node.age) ~ PoissonFirstEvent(lambda, age);
    walk(node.left, node.age, lambda);
    walk(node.right, node.age, lambda);
  }
  else
    # we waited until the present, but nothing happened.
   observe None ~ PoissonFirstEvent(lambda, age)
}

I suppose that this fails to share the exp(-lambda*(time-node.age)) between the two branches, but I'm not sure if that matters.

Just some thoughts -- I hope they are helpful.