Suppose we have that the distribution of a count $x$ is given by a binomial distribution: $x\sim \mathcal{B}(n,\theta)$. Further, suppose that the sample size $n\sim NB(\lambda,\kappa)$, where this represents a negative binomial distribution with mean $\lambda$ and variance $\lambda^2/\kappa$. I can write a Mathematica function that implements this parameterisation of the negative binomial as follows:
fCustomNB[λ_, κ_] := NegativeBinomialDistribution[κ, κ/(κ + λ)]
which has desired mean (variance also but omitting to keep question short):
FullSimplify@Mean[fCustomNB[λ, κ]]
λ
I can also use Mathematica to find the distribution of $x$ where I have marginalised $n$ out of the distribution (i.e. the distribution of $x|\theta, \lambda, \kappa$) as follows:
aDist = ParameterMixtureDistribution[
BinomialDistribution[n, θ],
n \[Distributed] fCustomNB[λ, κ]];
I can then evaluate the PDF of this using:
Assuming[{λ > 0, κ > 0, 0 < θ < 1, x >= 0}, PDF[aDist, x]]
which gives me this, somewhat heathen, function:
(1 - θ)^(-x + Ceiling[x]) θ^x κ^κ λ^ Ceiling[x] (κ + λ)^(-κ - Ceiling[x]) Binomial[Ceiling[x], x] Binomial[-1 + κ + Ceiling[x], -1 + κ] Hypergeometric2F1[ 1, κ + Ceiling[x], 1 - x + Ceiling[ x], (λ - θ λ)/(κ + λ)]
The part of this that is particularly difficult to translate into other languages (this distribution represents part of a larger model already coded up in Stan) is the Hypergeometric bit.
The bottom line is, can I simplify this expression in terms of more vanilla-like functions?
Edit: as finishing writing this question, I realised we can -- so I include my answer below in case it helps others!