I'm not sure if this is the best way to go about this, because I'm fairly new to Bayesian methods. I'm trying to model a process where the number of trials $n$ used in a binomial process is generated by a non-homogeneous Poisson process. I'm currently using PyMC to fit the Poisson part of this (example here, without the whole capping part), but I can't figure out how to integrate this with the Binomial. How can I take what I've generated with the Poisson process and use it in fitting a Binomial process? Or is there a better way to do this using similar methods?


Here's what I've tried:

import pymc as pm
import numpy as np

t = np.arange(5)

a = pm.Uniform(name='a', value=1., lower=0, upper=10)
b = pm.Uniform(name='b', value=1., lower=0, upper=10)

def linear(a=a, b=b):
    return a * t + b

N_A = pm.Poisson(mu=linear, name='N_A')
C = pm.Beta('C', 1, 1)
obs_A = pm.Binomial('obs_A', N_A, C, observed=True, value=np.array([0,1,4,3,7]))

mcmc = pm.MCMC([obs_A, C, N_A, a, b])

When I try to pull a sample, it throws the error

pymc.Node.ZeroProbability: Stochastic obs_A's value is outside its support, or it forbids its parents' current values.

I'm sure I'm formulating this incorrectly, but I'm not sure how.


It looks like the observed value is exceeding the $n$ parameter of the binomial distribution.

I suspect you'll continue running into this problem so long as your parameter space includes values for which $p(D) = 0$. For example, a binomial distribution with $n=3$ assigns zero probability to $k \gt 3$, so if your MCMC chain sampled a value from the Poisson distribution that corresponded to a higher observed value, the data would have zero probability.

You may want to change your model such that the data are possible under any combination of parameters. In particular, you may want to model your observed values using a distribution that has infinite support. (Naturally, this will change how parameters are interpreted.)

Depending on your data, you may be able to treat $n$ as observed. For example, say you were modeling the number of users who walk into a store and express interest in a given phone. In that case, you could observe both $n_t$, the count of persons walking into the store, and $k_t$, the number expressing interest. I can't speak to whether this is realistic in the application you're facing, but it would likely avoid the ZeroProbability error.

  • $\begingroup$ Ok, that makes sense. I'm not sure I understand what you mean by reformulate though. $\endgroup$ – user3704120 Jul 26 '15 at 1:04
  • $\begingroup$ @user3704120 That have been a poor choice of term on my part: what I mean is that you may have to choose different distributions. In particular, you might want to model the observed values using a distribution with infinite support, say negative binomial. (Naturally this will change how parameters are interpreted.) $\endgroup$ – Sean Easter Jul 26 '15 at 1:34
  • $\begingroup$ Ahh ok, now I get what you're saying. Thank you for the help! $\endgroup$ – user3704120 Jul 26 '15 at 14:20
  • $\begingroup$ @user3704120 Quite welcome! I've edited my answer to include the details about choice of distribution, and added another thought on whether $n$ is observed. $\endgroup$ – Sean Easter Jul 26 '15 at 14:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.