Math (and other) issues in the R = 0 case

[afmagee42]

We have some conflicting opinions in the code about $R = 0$.

• In nbbp_ep we error out if $R = 0$
• In .dnbbp_subcrit we have a special case for $R =0$ that I think we inherited from when this was cdcent/chainsR.

There is also a related bug in dnbbp where the if (!condition_on_extinction) should be if (!condition_on_extinction && any(is.infinite(x))). We don't need the mass of non-finite chains if there aren't any.

(Found when working to resolve #7)

Much ado about math

Mathematically things get a bit fuzzy at $R = 0$. From a mean/variance perspective, for the mean to be 0, the mass at all other values has to be 0, so the variance is 0, and you should have to have $k = \infty$, rather than a free $k$.

On the other hand, it seems sane enough to assert that at $R = 0$ with probability 1 the chain size is 1 (and with probability 0 any other size), making the extinction probability 1.

For the PMF, we use Blumberg and Lloyd-Smith (2013) Equation 9

At $R = 0$, this becomes (for chain size $j$)

$$\frac{\Gamma(kj + j - 1)}{\Gamma(kj)\Gamma(j + 1)} \frac{0^{j - 1}}{1^{kj + j - 1}}$$

which becomes, at $j = 1$

$$\frac{1}{\Gamma(2)} \frac{0^0}{1^1} = 0^0$$

while for $j >= 2$ we get 0 raised to a positive integer power, so the whole thing becomes 0, regardless of $k$.

R regards the log-PMF at 0 as NaN because it involves -Inf * 0.

Moving forward

I see two solutions:

• Assert that $R$ is strictly positive and remove the special case entirely.
• Assert that at $R = 0$ the chain size distribution is $Pr(j = 1) = 1$, $Pr(j >= 2) = 0$ explicitly in code without calling dnbinom.

[afmagee42]

Thoughts @swo (who may remember the story behind the special case), @SamuelBrand1, or @seabbs?

[SamuelBrand1]

When are you sampling "exactly" 0?

[afmagee42]

Great question, and I realize I left out context. This not something encountered in any MCMC, it's sort of encountered in maximum likelihood estimation for a bunch of chains of size 1. There it converges to very small $R$ and arbitrary $k$.

So, this is sort of an attempt to answer to "what is the MLE for all singleton chains?"

[swo]

I think the way to choose is what makes the likelihood surface smooth? Like, it seems that

$$ \lim_{R \to 0} P[j = 1] = 1 $$

so it would make conceptual sense and also wouldn't break anything mathematically to assert that $P[j=1]=1$ for $R=0$?