These functions are primarily designed for speed in simulation. Arguments are
not checked.
Suppose \(X_1 \mid G = g \sim \text{Poisson}(\mu g)\) and
\(X_2 \mid G = g \sim \text{Poisson}(r \mu g)\) where
\(G \sim \text{Gamma}(\theta, \theta^{-1})\) is the random item (subject) effect.
Then \(X_1, X_2 \sim \text{BNB}(\mu, r, \theta)\) is the joint distribution where
\(X_1\) and \(X_2\) are dependent (though conditionally independent),
\(X_1\) is the count outcome for sample 1 of the items (subjects),
\(X_2\) is the count outcome for sample 2 of the items (subjects),
\(\mu\) is the conditional mean of sample 1, \(r\) is the ratio of the
conditional means of sample 2 with respect to sample 1, and \(\theta\) is
the gamma distribution shape parameter which controls the dispersion and the
correlation between sample 1 and 2.
The likelihood is
$$
\begin{aligned}
L(r, \mu, \theta \mid X_1, X_2) = & \left( \frac{\theta^{\theta}}{\Gamma(\theta)} \right)^{n} \times \\
& \frac{\mu^{\sum{x_{1i}} + \sum{x_{2i}}}}{\prod_{i=1}^{n} x_{1i}!} \frac{r^{\sum{x_{2i}}}}{\prod_{i=1}^{n} x_{2i}!} \times \\
& \frac{\prod_{i = 1}^{n} \Gamma(x_{1i} + x_{2i} + \theta)}{(\mu + r \mu + \theta)^{\sum x_{1i} + x_{2i} + \theta}}
\end{aligned}
$$
and the parameter space is
\(\Theta = \left\{ (r, \mu, \theta) : r, \mu, \theta > 0 \right\}\).
The log-likelihood is
$$
\begin{aligned}
l(r, \mu, \theta) = \ & n \left[ \theta \ln \theta - \ln \Gamma(\theta) \right] + \\
& n (\bar{x}_1 + \bar{x}_2) \ln(\mu) + n \bar{x}_2 \ln r + \\
& \sum_{i=1}^{n}{\ln \Gamma(x_{1i} + x_{2i} + \theta)} - \\
& n (\bar{x}_1 + \bar{x}_2 + \theta) \ln(\mu + r\mu + \theta) - \\
& \sum_{i = 1}^{n}{\ln x_{1i}!} - \sum_{i = 1}^{n}{\ln x_{2i}!}
\end{aligned}
$$