In this notebook, I show how to train spatial capture-recapture (SCR) models in PyMC. SCR expands upon traditional capture-recapture by incorporating the location of the traps in the analysis. This matters because, typically, animals that live near a particular trap are more likely to be caught in it. In doing so, SCR links individual-level processes to the population-level, expanding the scientific scope of simple designs.
In this notebook, I train the simplest possible SCR model, SCR0 (Royle et al. 2013, chap. 5), where the goal is estimating the true population size \(N\). Similar to the other closed population notebooks, I do so using parameter-expanded data-augmentation (PX-DA). I also borrow the concept of the detection function from the distance sampling notebook.
As a motivating example, I use the ovenbird mist netting dataset provided by Murray Efford via the secr package in R. The design of the study is outlined in Efford, Dawson, and Robbins (2004) and Borchers and Efford (2008). In this dataset, ovenbirds were trapped in 44 mist nets over 8 to 10 consecutive days during the summers of 2005 to 2009.
%config InlineBackend.figure_format ='retina'import arviz as azimport matplotlib.pyplot as pltimport numpy as npimport pandas as pdimport pymc as pmimport pymc_extras as pmximport pytensor.tensor as ptimport seaborn as sns# only necessary on MacOS Sequoia# https://discourse.pymc.io/t/pytensor-fails-to-compile-model-after-upgrading-to-mac-os-15-4/16796/5import pytensorpytensor.config.cxx ='/usr/bin/clang++'# hyper parametersSEED =42RNG = np.random.default_rng(SEED)BUFFER =100M =200# plotting defaultsplt.style.use('fivethirtyeight')plt.rcParams['axes.facecolor'] ='white'plt.rcParams['figure.facecolor'] ='white'plt.rcParams['axes.spines.left'] =Falseplt.rcParams['axes.spines.right'] =Falseplt.rcParams['axes.spines.top'] =Falseplt.rcParams['axes.spines.bottom'] =Falsesns.set_palette("tab10")def invlogit(x):'''Inverse logit function'''return1/ (1+ np.exp(-x))def euclid_dist(X, S, library='np'):'''Pairwise euclidian distance between points in (M, 2) and (N, 2) arrays''' diff = X[np.newaxis, :, :] - S[:, np.newaxis, :]if library =='np':return np.sqrt(np.sum(diff **2, axis=-1))elif library =='pm':return pm.math.sqrt(pm.math.sum(diff **2, axis=-1))def half_normal(d, s, library='np'):'''Half normal detection function.'''if library =='np':return np.exp( - (d **2) / (2* s **2))elif library =='pm':return pm.math.exp( - (d **2) / (2* s **2))def exponential(d, s, library='np'):'''Negative exponential detection function.'''if library =='np':return np.exp(- d / s)elif library =='pm':return pm.math.exp(- d / s)# coordinates for each trapovenbird_trap = pd.read_csv('ovenbirdtrap.txt', delimiter=' ')trap_count, _ = ovenbird_trap.shape# information about each traptrap_x = ovenbird_trap.xtrap_y = ovenbird_trap.yX = ovenbird_trap[['x', 'y']].to_numpy()# define the state space around the trapsx_max = trap_x.max() + BUFFERy_max = trap_y.max() + BUFFERx_min = trap_x.min() - BUFFERy_min = trap_y.min() - BUFFER# scale for plottingscale = (y_max - y_min) / (x_max - x_min)# plot the trap locationsplot_width =2plot_height = plot_width * scalefig, ax = plt.subplots(figsize=(plot_width, plot_height))# plot the trapsax.scatter(trap_x, trap_y, marker='x', s=40, linewidth=1.5, color='C1')ax.set_ylim((y_min, y_max))ax.set_xlim((x_min, x_max))ax.annotate('44 nets\n30m apart', ha='center', xy=(55, -150), xycoords='data', color='black', xytext=(40, 30), textcoords='offset points', arrowprops=dict(arrowstyle="->", color='black', linewidth=1, connectionstyle="angle3,angleA=90,angleB=0"))# aestheticsax.set_aspect('equal')ax.set_title('Mist net locations')ax.grid(False)plt.show()
One difference between spatial and traditional (non-spatial) capture is the addition of the trap identifier in the capture history. Whereas a traditional capture history is [individual, occasion], a spatial capture history might be [individual, occasion, trap].
In the ovenbird example, I ignore the year dimension, pooling parameters across years, which allows for better estimation of the detection parameters. My hack for doing so is treating every band/year combination as a unique individual in a combined year capture history. This is easy to implement, creates an awkward interpretation of \(N\) (see below).
# ovenbird capture historyoven_ch = pd.read_csv('ovenbirdcapt.txt', delimiter=' ')# create a unique bird/year identifier for each individualoven_ch['ID'] = oven_ch.groupby(['Year','Band']).ngroup()occasion_count = oven_ch.Day.max()# merge the datasets, making sure that traps with no detections are includedovenbird = ( ovenbird_trap.merge(oven_ch[['ID', 'Net', 'Day']], how='left') [['ID', 'Day', 'Net', 'x', 'y']] .sort_values('ID') .reset_index(drop=True))ovenbird.head(10)
ID
Day
Net
x
y
0
0.0
1.0
2
-50.0
-255.0
1
1.0
9.0
20
-50.0
285.0
2
1.0
1.0
15
-50.0
135.0
3
2.0
6.0
17
-50.0
195.0
4
2.0
9.0
27
49.0
165.0
5
2.0
1.0
15
-50.0
135.0
6
2.0
1.0
14
-50.0
105.0
7
3.0
1.0
41
49.0
-255.0
8
3.0
3.0
39
49.0
-195.0
9
3.0
1.0
42
49.0
-285.0
Simulation
Before estimating the parameters, I perform a small simulation. The simulation starts with a core idea of SCR: the activity center. The activity center \(\mathbf{s}_i\) is the most likely place that you’d find an individual \(i\) over the course of the trapping study. In this case, I assume that activity centers are uniformly distributed across the sample space.
I compute the probability of detection for individual \(i\) at trap \(j\) as \(p_{i,j}=g_0 \exp(-d_{i,j}^2/2\sigma^2),\) where \(g_0\) is the probability of detecting an individual when it’s activity center is at the trap, \(d_{i,j}\) is the euclidean distance between the trap and the activity center, and \(\sigma\) is the detection range parameter.
# true population sizeN =150# simulate activity centersS_true = RNG.uniform((x_min, y_min), (x_max, y_max), (N, 2))# true distance between the trap and the activity centersd_true = euclid_dist(X, S_true)# detection parametersg0_true =0.025sigma_true =73# simulate the number of captures at each trap for each individualcapture_probability = g0_true * half_normal(d_true, sigma_true)sim_Y = RNG.binomial(occasion_count, capture_probability)# filter out undetected individualswas_detected = sim_Y.sum(axis=1) >0sim_Y_det = sim_Y[was_detected]n_detected =int(was_detected.sum())
Following Royle et al. (2013), Chapter 5, I first fit the version of the model where we assume that we know the true population size. In this case, I’m only estimating the detection parameters and the activity center locations.
# upper bound for the uniform prior on sigmaU_SIGMA =150with pm.Model() as known:# priors for the activity centers S = pm.Uniform('S', (x_min, y_min), (x_max, y_max), shape=(n_detected, 2))# priors for the detection parameters g0 = pm.Uniform('g0', 0, 1) sigma = pm.Uniform('sigma', 0, U_SIGMA)# probability of capture for each individual at each trap distance = euclid_dist(X, S, 'pm') p = pm.Deterministic('p', g0 * half_normal(distance, sigma))# likelihood pm.Binomial('y', p=p, n=occasion_count, observed=sim_Y_det )pm.model_to_graphviz(known)
Figure 2: Visual representation of the model where \(N\) is known.
NUTS sampler
Throughout these notebooks, I have been estimating parameters with the default NUTS sampler in PyMC. One of PyMC’s great strengths is that it allows you to select different samplers with the same model code. We saw one version of this in the occupancy notebook, where we had a Binary Gibbs Sampler for the discrete \(z_i\) state and a NUTS sampler for the continuous parameters. But PyMC goes one step further, allowing you to choose different NUTS samplers. These include numpyro, which is itself a package for MCMC that can run on GPUs. Additionally, there is the nutpie sampler that is ultimately written in Rust. In most cases, nutpie is much faster than the default PyMC sampler. That said, nutpie can slightly less flexible.
Since the SCR models in this notebook take a little longer to run, I will use the nutpie sampler throughout. The difference can be substantial. For instance, the unknown \(N\) model below takes 39 seconds with the default sampler, and 13 seconds with the nutpie sampler.
with known: known_idata = pm.sample(nuts_sampler='nutpie')
Figure 3: Trace plots for model where \(N\) is known. The true parameter values are shown by vertical and horizontal lines.
The trace plots show reasonable agreement between the true parameter values and the estimated values, although \(g_0\) appears to be overestimated.
Ovenbird density
Now, I estimate the density \(D\) for the ovenbird population. Like distance sampling, SCR can robustly estimate the density of the population, regardless of the size of the state space. The difference between the model above and this one is that we use PX-DA to estimate the inclusion probability \(\psi,\) and subsequently \(N.\) First, I convert the DataFrame to a (n_detected, n_traps) array of binomial counts.
def get_Y(ch):'''Get a (individual_count, trap_count) array of detections.'''# count the number of detections per individual per trap detection_counts = pd.crosstab(ch.ID, ch.Net, dropna=False)# remove the ghost nan individual detection_counts = detection_counts.loc[~detection_counts.index.isna()] Y = detection_counts.to_numpy()return YY = get_Y(ovenbird)detected_count, trap_count = Y.shape# augmented spatial capture histories with all zero historiesall_zero_history = np.zeros((M - detected_count, trap_count))Y_augmented = np.vstack((Y, all_zero_history))
As with the other closed models in this series, I will write the model in terms of the latent inclusion state \(z_i\). Interestingly, .
with pm.Model() as oven:# Priors# activity centers S = pm.Uniform('S', (x_min, y_min), (x_max, y_max), shape=(M, 2))# capture parameters g0 = pm.Uniform('g0', 0, 1) sigma = pm.Uniform('sigma', 0, U_SIGMA)# inclusion probability psi = pm.Beta('psi', 0.001, 1)# compute the capture probability distance = euclid_dist(X, S, 'pm') p = pm.Deterministic('p', g0 * half_normal(distance, sigma))# inclusion state z = pm.Bernoulli('z', psi, shape=M)# likelihood mu_y = z[:, None] * p pm.Binomial('y', p=mu_y, n=occasion_count, observed=Y_augmented)pm.model_to_graphviz(oven)
Figure 4: Visual representation of the ovenbird model using data augmentation.
Figure 5: Trace plots for the ovenbird model using data augmentation. Maximum likelihood estimates are shown by vertical and horizontal lines.
The estimates are quite close to the maximum likelihood estimates, which I estimated with the secr package in R.
Finally, I estimate density \(D\) using the results. As in the closed capture-recapture and distance sampling notebooks, I use the posterior samples of \(\psi\) and \(M\) to sample the posterior of \(N.\) This \(N,\) however, has an awkward interpretation because I pooled across the years by combining all the detection histories. To get around this, I compute the average annual abundance by dividing by the total number of years in the sample. Then, I divide by the area of the state space.
N_samps = sim_N(oven_idata, detected_count, occasion_count)# kludgy way of calculating avergage abundanceyear_count =5average_annual_abundance = N_samps // year_count# area of the state space in terms of hectaresha =100*100mask_area = (x_max - x_min) * (y_max - y_min) / ha# densityD_samples = average_annual_abundance / mask_areaD_mle =1.262946fig, ax = plt.subplots(figsize=(4,4))ax.hist(D_samples, edgecolor='white', bins=13)ax.axvline(D_mle, linestyle='--',color='C1')ax.set_xlabel('Ovenbirds per hectare')ax.set_ylabel('Number of samples')ax.text(1.4, 800, rf'$\hat{{D}}$={D_samples.mean():.2f}', va='bottom', ha='left')plt.show()
Figure 6: Posterior distribution of the density \(D\) of ovenbirds. The maximum likelihood estimate is shown by the dotted red line.
Sometimes, the location of the activity centers is of interest. Below, I plot the posterior median for the activity centers for the detected individuals.
Figure 9: Posterior distribution for the detection function. The line represents the posterior mean while the shaded area is the 96% interval.
%load_ext watermark%watermark -n -u -v -iv -w
The watermark extension is already loaded. To reload it, use:
%reload_ext watermark
Last updated: Tue Oct 14 2025
Python implementation: CPython
Python version : 3.13.3
IPython version : 9.2.0
pytensor : 2.31.7
pandas : 2.2.3
matplotlib : 3.10.1
arviz : 0.21.0
seaborn : 0.13.2
numpy : 2.2.5
pymc_extras: 0.5.0
pymc : 5.25.1
Watermark: 2.5.0
References
Borchers, David L, and MG2432407 Efford. 2008. “Spatially Explicit Maximum Likelihood Methods for Capture–Recapture Studies.”Biometrics 64 (2): 377–85.
Efford, Murray G, Deanna K Dawson, and Chandler S Robbins. 2004. “DENSITY: Software for Analysing Capture-Recapture Data from Passive Detector Arrays.”Animal Biodiversity and Conservation 27 (1): 217–28.
Royle, J Andrew, Richard B Chandler, Rahel Sollmann, and Beth Gardner. 2013. Spatial Capture-Recapture. Academic press.