Building an abudance and detection probability model in PyMC3

Questions
1

Hello,

I’m trying to build an abundance and detection probability model in PyMC3 and have gotten stuck in the PyMC3 syntax - I hope someone can help me out. Thanks in advance!

Here’s a toy problem to demonstrate the issue:

I want to model the number of birds in my yard on any given day. I’ve counted the birds in my yard every day for 100 days, but I’m not a perfect birder, meaning that there could be birds present that I didn’t see. There are two parts to the model - the first is a Poisson model for the counts of birds present (abundance, b). The second is a Binomial model for the detection probability, which is the probability that I see a bird, given that it’s present in my yard (detection probability, p). Using my counts of birds over 100 days, I want to simultaneously estimate the abundance (b) and the detection probability (p).

obsbirds_{d} \sim Binomial(birds_{d},p)
birds_{d} \sim Poisson(b)

import pymc3 as pm
import scipy.stats as stats

#generate data 

b = 1 #Poisson rate parameter (average number of birds per day in my yard)
p = 0.3 # detection probability (probability that I see a bird, given that it's present)
d = 100 # number of days that I watched birds

birds = stats.poisson.rvs(b, size = d)
obsbirds = stats.binom.rvs(birds, p)

I now want to build a model in PyMC3 that estimates b and p simulataneously. Here’s my attempt, but the chains are failing and I can’t figure out what’s wrong.

# instantiate pymc3 model
mymodel = pm.Model()

# set up
with mymodel:

    ##priors 
    
    #for poisson part of model

    bfit = pm.HalfNormal('b', sigma = 10, shape = 1) #prior for b (must be positive)
    
    #for binomial part of model
    
    pfit = pm.Uniform('p', lower = 0, upper = 1, shape = 1) #prior for p (must be between 0 and 1)
    
    ##likelihoods
    
    birdsfit = pm.Poisson('birds', mu = bfit, shape = len(obsbirds))
    
    
    countsobs = pm.Binomial('obscounts',
                            n = birdsfit, 
                            p = pfit,
                            observed=obsbirds)

I’d appreciate help with the PyMC3 syntax, or pointers in the right direction. Thanks!

1 Like
2

There’s some trick with shape to get these likelihoods to work with n ~ (d,), p ~ 1, observed ~ (d,). But weirdly this doesn’t work:

 countsobs = pm.Binomial('obscounts',
                            n = birdsfit, 
                            p = pfit,
                            observed=obsbirds,
                            shape=len(obsbirds))

nor does direct indexing (?!?!)

''
for i in range(len(obsbirds)):
    co_i = pm.Binomial('co_%d' % i,
                      n=birdsfit[i],
                      p=pfit,
                      observed=obsbirds[i])

both give the same shape mismatch error.

1 Like
3

Hi :slight_smile:

I’ve just ran your code snippet and I was able to sample, here the traceplots (the model did not converge)

trace
trace872×440 182 KB

Could you specify which kind of error are you getting?

1 Like
4

Interesting, I don’t sample under python 3.8.8 and pymc3 3.11.0 (Ubuntu 20.04)

5

Thanks for trying this out! Interesting that it’s running for @vb690 .

With Python 3.7.0, pymc3 3.8, and Ubuntu 18.04, I’m getting this error:

image
image1360×481 61 KB

6

Hello

I run
Python Version 3.8.8
PyMC3 Version 3.11.2 installed through conda 4.9.2

On Windows 10 (yes for one time Windows OS is not the problem).

@lvluna Maybe you could try to update both Python and PyMC3 (@chartl only PyMC3), I do not know if it is related but in the release notes for PyMC3 3.11.1 they have Automatic imputations now also work with ndarray data, not just pd.Series or pd.DataFrame (see#4439).

Let me know if it changes something.

7

In looking at the traceplots that you get, I’ve realized that there is something wrong in the way that I’ve specified the model, not just in whether the sampler is running. I’m going to have to think about this, but will update when I figure it out. Thanks!

1 Like