In [1]:
%pylab inline
import ABCPRC as prc
import seaborn as sns
import scipy.stats as stats
First create the Stochastic SIS model. Briefly, for $n$ individuals, each susceptible individual has a rate of becoming infected $\beta I/n$ and each infected individual can recover at rate $\gamma$.
Here we use a tau-step algorithm to generate a realisation of the SIS process, where for each time-step a Poisson number of infection and recovery events occur according to a Poisson distribution. A small fixed time-step eps is also defined.
The function createData generates a realisation of the stochastic process. The function ibm is what is fed into the ABC method. It does a transformation of the parameters from the infection time and $R_0$ to $\beta$ and $\gamma$. Importantly the output of the ibm function is a two-dimensional time-series representing the $I_{t+1} | I_t$ pairs from the realisation.
It also does a small amount of error handling for when the infected population dies out.
In [34]:
def createData(beta,gamma,n=100,T=100,I0=1):
S,I = np.zeros(T),np.zeros(T)
S[0] = n-I0
I[0] = I0
eps=0.1
for i in np.arange(1,T):
if beta*S[i-1]*I[i-1]/n<0: print(beta*S[i-1]*I[i-1]/n)
infections = np.random.poisson(lam=eps*beta*S[i-1]*I[i-1]/n)
recoveries = np.random.poisson(lam=eps*gamma*I[i-1])
S[i] = np.max((0,S[i-1] - infections + recoveries))
I[i] = np.max((0,I[i-1] + infections - recoveries))
return S,I
def ibm(*ps):
n=100
gamma = 1./ps[1]
beta = ps[0]*gamma
S,I = createData(beta,gamma,n=n,T=100,I0=20)
if I[-1]==0:
return np.vstack((0*I[:-1],0*I[1:])).T
else:
return np.vstack((I[:-1],I[1:])).T
In [41]:
r0 = 1.5
rec_time = 1.
gamma = 1./rec_time
beta = r0*gamma
xs = ibm(r0,gamma)
Check the shape of the output. $T=100$, so the number of pairs should be $T-1$ and the second dimension should be $2$
In [42]:
print(xs.shape)
In [43]:
plt.plot(xs[:,0]);
In [32]:
plt.hist(stats.norm(loc=5.,scale=1.).rvs(2000));
In [78]:
m = prc.ABC()
We'll begin by misspecifying priors by having a uniform prior for $R_0$ and the revovery time that is a long way away from the true values. An error should raise that the prior is misspecified
In [72]:
wrong_priors = [stats.uniform(loc=0.,scale=0.1).rvs,stats.uniform(loc=5.,scale=0.1).rvs]
m.setup(modelFunc=ibm,xs=xs,priors=wrong_priors,method='Adaptive',toln=12)
m.run(10)
Let's try that again with properly specified priors that include the true value.
Here we define the method by calling the ABC class from the prc library. Priors are specified using the scipy stats package. In general priors are expected to be a list of functions where each function call returns a sample from the given parameter. The length of the prior list is used to determine the size of the inputs to the model.
The method is set-up by specifying the model modelFunc, the data xs, the priors priors, the method and the number of tolerances to use in the method toln.
The method is then ran using the method run for 100 particles.
In [84]:
m = prc.ABC()
priors = [stats.uniform(loc=0.,scale=2.).rvs,stats.uniform(loc=0.5,scale=2.).rvs]
m.setup(modelFunc=ibm,xs=xs,priors=priors,method='Adaptive',toln=10)
m.run(100)
Print out the summary of the parameters (median and 95% lower and upper percentiles)
In [85]:
summary = m.fitSummary()
Plot the distribution of the particles
In [86]:
m.trace(plot=True)