gistfile1.py

#observed data
popsize = 10000.0 
infected_data = np.array([  7, 10, 11, 11, 13, 14, 14, 13, 14, 16, 18, 18, 20, 21, 22, 23, 23, 24, 25, 23])
T = len(infected_data)

I0 = 7
R0 = 5
 

# priors
xbeta = pymc.Uniform('xbeta', 0., 1., value=0.5)
xgamma = pymc.Uniform('xgamma', 0., 1., value=0.5)
S0 = pymc.Uniform('S0', 0, popsize, value=36)

def crop(x, minval, maxval):
    return min(max(x, minval), maxval)

# deterministic compartmental model
def simulate(S0, beta, gamma, timespan=T):
    S = np.zeros(timespan)
    I = np.zeros(timespan)
    R = np.zeros(timespan)
    S[0] = S0
    I[0] = I0
    R[0] = R0
    for i in range(1,timespan):
        newly_infected_count = crop(beta*S[i-1]*I[i-1]/popsize, 0., S[i-1])
        recovered_count = crop(gamma*I[i-1], 0., I[i-1])
        
        S[i] = S[i-1] - newly_infected_count
        I[i] = crop(I[i-1] - recovered_count + newly_infected_count, 0., popsize)
        R[i] = crop(R[i-1] + recovered_count, 0., popsize)    
    return S, I, R

@pymc.deterministic
def SI(S0=S0, xbeta=xbeta, xgamma=xgamma):
    s, i, r = simulate(S0, xbeta, xgamma)
    return i
    
infected = pymc.Poisson('infected', mu=SI, value=infected_data, observed=True)


#run the model
np.random.seed(1701)
model=pymc.Model([xbeta, xgamma, S0, SI, infected])
map_ = pymc.MAP( model )
map_.fit()
mc = pymc.MCMC(model)
mc.sample(iter=500, burn=100)

# weird result
assert(xbeta.value > 1)