0

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.

Community
  • 1
  • 1
user3704120
  • 141
  • 1
  • 8

1 Answers1

0

The error appears because the automatically generated initial value of N_A (the number of trials) is less than some of the observed values specified in pm.Binomial, i.e. the probability of observing such a number is zero, and the initial value has to be rejected. One solution would be to supply an acceptable initial value for N_A explicitly:

x = np.array([0,1,4,3,7])
N_A = pm.Poisson(mu=linear, name='N_A', value=x)

For reference: https://github.com/pymc-devs/pymc/issues/12:

This is a model mis-specification issue, not a bug. It means that the initial parameter value for the likelhood are incompatible with the data -- you need to specify compatible initial values.

peter
  • 480
  • 3
  • 7