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) @pm.deterministic 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]) mcmc.sample(10000,5000)
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.