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

I asked this question in stackoverflow but closed it after 4 days since I didn't get any answer. I am using this example as a basis for my work: I am trying to use pymc MCMC to fit the parameters of a deterministic model. THe model has 8 equations and 7 parameters, and describes an outbreak of an infectious disease transmitted through vectors (mosquitoes). The data I have are the number of new cases in any given day, but only a fraction of the infections become cases, and out of those, only a fraction of them are reported. So, for the days that I have no cases reported, I replaced the 0 entries with "None" so that PYMC knows that these are missing values.

So, I have an ODE system that models the transmission and dynamics of the outbreak with 7 unknown parameters, and 2 additional parameters for the fraction of cases ($symFraction$) and the fraction of reported cases ($repRate$).

The idea is that the ODE gives me the count of the number of new infections, and

data ~ NegativeBinomial

with $mean = repRate*symFraction*new infections$.

I gave priors to some of the parameters of interest:

#priors:
betah = pc.Uniform('betah', 0.0, 1e6, value= 0.001) #force of infection for humans
betamos = pc.Uniform('betamos', 0.0, 1e6, value= 0.001) #force of infection for mosquitoes
mos2humanRatio = pc.Uniform('mos2humanRatio', 0, 1e6, value=10.0) #mosquito density
symfraction = pc.Uniform('symfraction', 0.0, 1.0, value=0.2) #proportion of infections that become cases
repRate = pc.Uniform('repRate', 0.0, 1.0, value=0.1) #proportion of cases reported
disper_par = pc.Uniform('disper_par', 0, 1e6, value= 1)


and I defined my likelihood as:

A = pc.NegativeBinomial('A', mu=incidence, alpha=disper_par, value=masked_values, observed=True)


My problem is that I keep getting this error for my likelihood:

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


I have read all over the place to try to understand this, and I just cannot fix it. So I guess I have 3 questions, one statistical, the other 2 for pymc experts:

1. Is this the correct way of setting up the MCMC?
2. Am I handling the missing data correctly?
3. Why I keep getting this error? I just don't understand why the likelihood is outside its support... or why I would not allow the values of the parents...

I don't quite understand if the problem resides with the missing data or with one of the priors... I am completely new at pymc and MCMC, so any help would be highly appreciated!

The following is a mwe:

from __future__ import division
import pymc as pc
import numpy as np
from matplotlib import pyplot as plt
from scipy.integrate import odeint

def my_eqs(y,t,params):
#parameters:
[betah, betamos, bmos, delta,  gammamos, mu, N, psi] = params
S_h, E_h, Inf_h, R_h, C_h = y[0], y[1], y[2], y[3], y[4]
S_mos, E_mos, Inf_mos = y[5], y[6], y[7]
#human equations:
dS_h = -betah*(Inf_mos)*S_h
dE_h = betah*S_h*Inf_mos - delta*E_h
dInf_h = delta*E_h - psi*Inf_h
dR_h = psi*Inf_h
dC_h = delta*E_h

#mosquito equations:
dS_mos = bmos -betamos*(Inf_h/N)*S_mos - mu*S_mos
dE_mos = betamos*(Inf_h/N)*S_mos - mu*E_mos -gammamos*E_mos
dInf_mos = gammamos*E_mos -mu*Inf_mos

dydt = [dS_h, dE_h, dInf_h, dR_h, dC_h, dS_mos, dE_mos, dInf_mos]
return dydt

#create some fake data:
mydata = [0.0, 0.0, 0.0, 3.0, 8.0, 7.0, 1.0, 5.0, 3.0, 1.0, 0.0, 6.0, 9.0, 2.0, 0.0, 3.0, 8.0, 8.0, 6.0, 0.0, 2.0, 4.0,
4.0, 9.0, 8.0, 6.0, 8.0, 0.0, 2.0, 1.0, 0.0, 7.0, 5.0, 9.0, 1.0, 4.0, 6.0, 7.0, 4.0, 10.0, 4.0, 2.0, 12.0, 6.0, 4.0,
5.0, 9.0, 12.0, 10.0, 10.0, 12.0, 8.0, 15.0, 17.0, 15.0, 30.0, 32.0, 39.0, 44.0, 51.0, 57.0, 61.0, 31.0, 32.0, 25.0,
39.0, 45.0, 54.0, 65.0, 56.0, 60.0, 57.0, 48.0, 47.0, 66.0, 48.0, 39.0, 48.0, 32.0, 39.0, 45.0, 41.0, 59.0, 38.0, 41.0,
41.0, 48.0, 40.0, 35.0, 38.0, 27.0, 29.0, 20.0, 14.0, 17.0, 17.0, 16.0, 11.0, 10.0, 11.0, 22.0, 21.0, 6.0, 12.0, 9.0,
14.0, 4.0, 7.0, 15.0, 14.0, 13.0, 6.0, 12.0, 49.0, 22.0, 9.0, 6.0, 8.0, 0.0, 5.0, 12.0, 5.0, 10.0, 8.0, 11.0, 15.0, 5.0,
9.0, 6.0, 3.0, 3.0, 2.0, 2.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

#replace the days where mydata=0 with None to handle missing data:
cases = np.array(mydata, dtype=object)
cases[1:][cases[1:] == 0] = None

# plt.plot(cases,'o')

###### MCMC setup
#deterministic parameters:
tspan = range(0,len(mydata),1)
delta = 1/4.5
gammamos = 0.1
mymu = 1/15.0
psi = 1/5.5
N = 100000

#priors:
betah = pc.Uniform('betah', 0.0, 1e6, value= 0.001)
betamos = pc.Uniform('betamos', 0.0, 1e6, value= 0.001)
mos2humanRatio = pc.Uniform('mos2humanRatio', 0, 1e6, value=10.0)
symfraction = pc.Uniform('symfraction', 0.0, 1.0, value=0.2) #proportion of infections that become symptomatic
repRate = pc.Uniform('repRate', 0.0, 1.0, value=0.1) #proportion of cases reported
disper_par = pc.Uniform('disper_par', 0, 1e6, value= 1)

@pc.deterministic
def bmos(mos2humanRatio=mos2humanRatio):
out = mos2humanRatio*(N)*mymu
return out

####deterministic model
@pc.deterministic
def mymodel(repRate=repRate,symfraction=symfraction, betah=betah,
betamos=betamos,mos2humanRatio=mos2humanRatio,bmos=bmos):
initCond = [N-1, 0, 1, 0, 0, mos2humanRatio*N, 0, 0]
params = [betah, betamos, bmos, delta,  gammamos, mymu, N, psi]
soln = odeint(my_eqs, initCond, tspan, args = (params, ))
newcases = np.zeros(np.size(soln[:, 4]))
for t in range(1,len(tspan)):
newcases[t] = (soln[t, 4] - soln[t-1, 4])
reported_cases = repRate*symfraction*newcases
# plt.plot(reported_cases, 'r')
return reported_cases

incidence = pc.Lambda('incidence', lambda mymodel=mymodel: mymodel)

A = pc.NegativeBinomial('A', mu=incidence, alpha=disper_par, value=masked_values, observed=True)

• Welcome to Cross Validated. It seems you have a statistical question somewhere in there buried under a ton of code. Questions that are solely about programming are off-topic for this site and may be removed. Can you edit your question to emphasize your statistical question? For example, explain a bit more about the model you are trying to estimate and give more context about the data that you have? I.e., explaining what the outcome(s) and explanatory variables are – Marquis de Carabas May 4 '16 at 21:13
• @marquisdecarabas: thanks for your comment. I edited the question and hopefully now it is more reflective of my question(s)... – Laura May 4 '16 at 21:56