N-mixtures Transition Probability Matrix

The likelihood function for the open-population N-mixtures model from Dail and Madsen (2011) is given in Eq. (5). In this model there are four parameters: probability of detection p, initial site abundance λ, mean population growth rate γ, and survival probability ω. There are also two study specific constants: the number of sampling sites R, and the number of sampling occasions T. The total population at site i and time t are given by the latent variables Nit, and the observed population counts are denoted nit.

The function Pa,b shown in Eq. (6) calculates the transition probability for moving from a population of size a to a population of size b. To calculate the likelihood equation requires calculating Pa,b at most RT(K+1)(T-1) times (when each nit=0), where K is the upper bound on the summations and thus also the upper bound on the population size. The complexity of computing Pa,b is thus the main bottleneck for computing the likelihood function. Pa,b is a convolution of two discrete distributions and can thus be calculated efficiently using the FFT convolution.

Here Bin and Pois denote the binomial and Poisson distribution functions, respectively. We define the transition probability matrix MK to be the matrix of Pa,b values, where a and b vary from 0 to K.

Equation (7) illustrates the relationships between Pa,b, MK, convolution, and FFT; let xa,c=Bin(c;a,ω) so that x_a={xa,0,xa,1,…,xa,min{a,b},0,0,…,0} (right padded with zeroes until x_a has K+1 elements), and let yb,c=Pois(b-c;γ) so that y_={y0,c,y1,c,…,yK,c}, then (x_a∗y_)b=∑c=0min{a,b}xa,c·yb,c=Pa,b, and (x_a∗y_)={(x_a∗y_)0,(x_a∗y_)1,…,(x_a∗y_)K}.

Computation of MK is a large portion of the computational cost of calculating the likelihood function in open-population N-mixture models. In Fig. 1 we show a plot of the median computation times for calculating MK for increasing values of K using both the manual convolution and the FFT convolution techniques. The values plotted are the median computing times over 100 runs for each value of K. Each run involved generating a random value for ω drawn from a Beta (a=10, b=10) distribution and a random value for γ drawn from a χK2 distribution. MK was then calculated using the same ω and γ for both the manual convolution method and the FFT convolution method, and computation times were calculated using the R package microbenchmark (Mersmann 2021). The error bars shown in Fig. 1 indicate minimum and maximum observed computing times. We also note that for both methods, the upper and lower quantiles are indistinguishable from the median computing times on the log plot.

Plot of log median value of computation time measured in milliseconds (out of 100 runs with randomised γ and ω) versus K for computation of the transition probability matrix MK using manual (grey) and the fast Fourier transform (black) convolution methods. Error bars indicate minimum and maximum computing time over the 100 runs per K. We considered K= 1, 10, 50, 100, 250, 500, and 1000

Often, time-varying covariates for the population dynamics parameters γ and ω are incorporated into N-mixture models. When this is done, the matrix MK is calculated once for each time point t∈{2,…,T}. Thus the computing time savings when using FFT are increased by a factor of up to T-1 when time covariates are considered.

In Sect. 3.1 we use a simulation study to illustrate the improved computational efficiency of using the FFT method for calculating MK when fitting N-mixture models. In Sect. 4.1 we apply the FFT method to improve computing efficiency for an Ancient Murrelet chick counts application.