<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Stochastic SIS ABC fitting example
Step1:  Create individual-based model
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.
Step2:  Generate some fake data
Use the ibm function to create data for a specified $R_0$ and recovery time.
Step3:  Check the shape of the output. $T=100$, so the number of pairs should be $T-1$ and the second dimension should be $2$
Step4:  Set-up model
Step5:  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
Step6:  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.
Step7:  Print out the summary of the parameters (median and 95% lower and upper percentiles)
Step8:  Plot the distribution of the particles
<ASSISTANT_TASK:>
Python Code:
%pylab inline
import ABCPRC as prc
import seaborn as sns
import scipy.stats as stats
Explanation: Stochastic SIS ABC fitting example
End of explanation
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
Explanation: Create individual-based model
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.
End of explanation
r0 = 1.5
rec_time = 1.
gamma = 1./rec_time
beta = r0*gamma 
xs = ibm(r0,gamma)
Explanation: Generate some fake data
Use the ibm function to create data for a specified $R_0$ and recovery time.
End of explanation
print(xs.shape)
plt.plot(xs[:,0]);
plt.hist(stats.norm(loc=5.,scale=1.).rvs(2000));
Explanation: Check the shape of the output. $T=100$, so the number of pairs should be $T-1$ and the second dimension should be $2$
End of explanation
m = prc.ABC()
Explanation: Set-up model
End of explanation
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) 
Explanation: 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
End of explanation
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) 
Explanation: 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.
End of explanation
summary = m.fitSummary()
Explanation: Print out the summary of the parameters (median and 95% lower and upper percentiles)
End of explanation
m.trace(plot=True)
Explanation: Plot the distribution of the particles
End of explanation
<END_TASK>