0
$\begingroup$

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')

masked_values = np.ma.masked_equal(cases, value=None)


###### 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)
$\endgroup$
  • 1
    $\begingroup$ 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 $\endgroup$ – Marquis de Carabas May 4 '16 at 21:13
  • $\begingroup$ @marquisdecarabas: thanks for your comment. I edited the question and hopefully now it is more reflective of my question(s)... $\endgroup$ – Laura May 4 '16 at 21:56
2
$\begingroup$

It looks that the observed values exceeds the dominion of the distribution. When a random variable is defined as a function of another random variable, PyMC checks that no value of the parent distribution leads to an impossible value for the child distribution. Please look here for a more comprehensive explanation on a similar subject.

$\endgroup$
  • $\begingroup$ I wonder if you mean a slightly different word when you talk about the "dominion" of the distribution - did you intend "domain"? $\endgroup$ – Silverfish Jul 27 '16 at 18:27

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.