Hello,

In my model, I have a term of the following form: \frac{1-e^x}{1-k e^x}, where k>0 and x can range from -inf to inf.

In my case, such a term appears in several closed form expressions describing the simple birth and death process (a.k.a. Kendall process, a Markov birth death process). For example, the probability of the process going extinct by time t is P_{extinct}[t;\lambda,\mu]=\frac{\mu}{\lambda} \left( \frac{e^{(\lambda-\mu)t}-1}{e^{(\lambda-\mu)t}-\frac{\mu}{\lambda}} \right).

I am running into issues with numerical precision when the x in \frac{1-e^x}{1-k e^x} approaches -inf. (In my case, x=(\lambda-\mu)t and t is my data in minutes. Some values of t reach ~80 minutes, which makes x take on values like -60.) The term approaches 1, but I cannot have the term equal 1 exactly.

In one case I’ve debugged, I can’t have the term equal to 1 because I am using it as a probability `p_extinct`

in a simplex: `simplex theta = [prob_extinct, prob_event2, prob_event3]`

. The components of `theta`

must sum up to 1. When I calculate `prob_extinct`

=1 and `prob_event2`

= 1e-80, for example, then attempt `prob_event3 = 1 - prob_extinct - prob_event2`

, then I get `prob_event3`

< 0, which leads to a NaN if I work with log probabilities. Alternatively, if I compute `prob_event3`

separately, Stan throws an error because the sum of the terms of `theta`

is >1. My Stan model fails to initialize when run because the initial 100 parameter guesses (e.g. for \lambda and \mu) always lead to `theta`

containing NaN’s.

Computing \log(\frac{1-e^x}{1-k e^x}) as `log1m_exp(x) - log1m_exp( x + log(k) ) `

helps with precision of the numerator and denominator separately, but doesn’t get around the fact that I’m subtracting two similar numbers from each other. Is there another clever way to compute \frac{1-e^x}{1-k e^x} using the composed functions available in Stan?

Or is the problem elsewhere, perhaps in that I am not using simplex variables in their intended usage? Should I have an if statement like `if(p_extinct == 1){ return(1 - machineeps()) }`

to make sure I never return exactly 1 for `p_extinct`

?

In my model, the simplex is used like this: `mycountdata ~ multinomial( theta )`

.