{'fraction_non_alphanumeric': 0.030729132444094004, 'fraction_numerical': 0.02122347066167291, 'mean_word_length': 4.371529411764706, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 5, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stochastic Epidemic Modelling

Inferring how an epidemic will progress and what actions to take when presented with limited information is of critical importance for epidemiologists and health professionals. In real world settings, epidemiology data can be scarce or subject to reporting errors. In this project there are different epidemic scenarios simulated and, using hidden Markov Chains, it is attempted to mimic the imperfect data an epidemiologist will encounter. Furthermore, different kinds of compartmental models are modelled using the particle Markov Chain Monte Carlo algorithm with a variation of the adaptive Metropolis-Hastings algorithm to estimate the posterior density of the parameters underlying the models. Moreover, the sensitivity of these algorithms is investigated when subjected with changes in the dataset. This is accomplished by limiting the information provided, while using an adaptive approach on the posterior covariance of the parameters.


Introduction
In an epidemic setting it is common to be interested in the parameters underlying the disease, in order to be able to make health policy decisions on how to contain the epidemic or create a timeline of how it is expected to evolve. Compartmental models are a relatively simple and quick way to get some estimates for the disease's characteristics on a population level and make informed decisions on how to react, which were also used in the beginning of the COVID-19 pandemic ( [1], [25]). In this project the use of the Particle Markov Chain Monte Carlo (PMCMC) algorithm is explored, in modelling these stochastic epidemic models. The use of PMCMC in epidemic modelling has been investigated before, but not in the same manner as here. Dureau and others [8] researched the use of PMCMC in the SEIR model, with some variations, and time-varying parameters. Rasmussen and others [21] created a framework relying on PMCMC to combine compartmental models and genetics data to obtain the underlying disease dynamics and more recently Endo and others [10] presented a paper that acts as a tutorial for anyone who wants to apply the PMCMC algorithm to model epidemics with compartmental models. In Section 3, the methodology is presented -the theory behind the compartmental models, specifically the SIR model, as well as the algorithms used to simulate stochastic epidemic data and estimate their parameter's posterior densities. In Section 4, the results of fitting different scenarios using Approximate Bayesian Computation (ABC) and PMCMC are presented, namely results from modelling SIR and SEIR compartmental models with hidden Markov Chains to simulate real world scenarios. Two real world scenarios are looked into, noisy data and underreported data which are simulated with the Normal and the Binomial distributions respectively. After that, some experiments are reported on whether it is still possible to obtain Markov Chains that provide us with good posterior samples, while limiting the information in the dataset and reducing the amount of data the model trains on. It is attempted to limit the simulated data by increasing the noise in the noisy data and decreasing the probability of observation of cases in the under-reported data. Moreover, also investigated is what happens when trying to infer the parameters underlying those models earlier in the disease's progression, as it is imperative for healthcare professionals to be able to understand the disease's characteristics as early as possible. Finally, in the last section, Section 5, the possible extensions of this work are discussed as well as ideas worked on in this project that did not provide any meaningful results. Additionally, mentioned are some of the limitations of this algorithm in an epidemiology setting and some of the difficulties encountered.

Epidemiological Compartmental Models
When dealing with epidemics that occur within a limited time-frame, so that outside forces affecting the population (births and deaths) are negligible, a common and simple model to use is the SIR without demographics by Kermack WO and McKendrick AG [16]. This simple model can be split in 3 compartments: S which stands for Susceptible, the population that can contract the disease, I which stands for Infected, the population that can infect others, and R which stands for Removed, the population that has either recovered or died of the disease and ,thus, can no longer contract it. In this simple setting we assume that all contacts between any compartment are homogeneous and that people who belong in the same compartments adhere to the same disease dynamics. Hence we say that at time t, S(t), I(t), R(t) define the number of people who belong in each compartment. Since our model is a closed system we have that where N is the total population, which is constant across time. The system is defined by a set of three equations: Where β represents the rate at which people leave the Susceptible class to enter the Infected class and γ the rate at which people leave the Infected class to enter the Removed class. A way to interpret β and γ in the epidemics setting is that β is the rate of contact of a susceptible person times the probability of infection and 1/γ is the average time it takes an infected person to recover or die from the disease.
The last line in equation (1) is redundant, since we can simply find R(t) at each time-step from the other two equations using R(t) = N − I(t) − S(t). Hence the system can be described by: We say that the disease has run its course when at t we have S(t ) = 0, meaning that everybody susceptible to the epidemic has been actually infected, or when I(t ) = 0, meaning there are no more infected individuals. By looking at the equations above we can see that I (t) > 0 if and only if βS(t) − γ > 0. Since we need I (0) > 0 in order to have an epidemic we need S(0) > γ/β. We define R 0 := βN/γ, when I(0) is very small compared to S(0), as the basic reproduction number which is critical to whether the disease is going to progress to an epidemic or it is going to die out. If R 0 < 1 the disease will finish progressing with people still not infected, while if R 0 > 1 the disease will finish progressing when everybody has been infected.
There are many different compartmental models used depending on the nature of the disease. Another model I will investigate is the SEIR model. In the SEIR there is an extra compartment, E, namely the Exposed population. If a person is in the exposed compartment they have contracted the disease, but they are not yet infectious, meaning they cannot infect others. Since our model remains a closed system, births and deaths of natural causes are insignificant in the time it takes the epidemic to progress, we have S(t) + E(t) + I(t) + R(t) = N for all t, The SEIR model is defined by the set of the following equations: In the SEIR model we have an extra parameter, α, which is the rate of leaving the exposed compartment to enter the infectious one. It can be interpreted as 1/α signifying the average period a person spends in the exposed class, or the latent period. One thing to be noted about the SEIR model is that to enter the infectious compartment, you first need to enter the exposed one, so you cannot be infectious immediately after you contact the disease.
Having defined the deterministic SIR model we can now define the stochastic SIR model which progresses in continuous time [13]. The system can now be described by a Markov jump process in continuous time where there are two types of jumps, (S, I) → (S − 1, I + 1), when a person gets infected, and (S, I) → (S, I − 1), when a person recovers or dies. These events occur at a categorical probability with the times in between each jump being identically and independently exponentially distributed. Thus, in a time interval [t, t+∆t] we have the following probabilities of a jump happening that define the stochastic SIR model: The process (S t , I t ) for t ≥ 0 is a Poisson process where the stochastic rate λ(t) depends on time and the state of the process at time t. The jump (S, I) → (S − 1, I + 1) occurs at a stochastic rate β StIt N and the jump (S, I) → (S, I − 1) occurs at a stochastic rate γI t . We can clearly see that the probabilities for a jump to occur at a time interval [t, t + ∆t] are conditional at the state of the system at time t.

Gillespie Algorithm
To simulate a stochastic SIR model I will be using the Gillespie Algorithm. The algorithm, created by Jacob L. Doob [5] and popularised by Daniel T. Gillespie in 1976 [11], is used to simulate stochastic trajectories of a system where the rates of the random variables are known. Before we define the algorithm, we need to define the propensity functions of a stochastic system. These are the stochastic rates of the jumps in the system that were defined previously.
The direct formulation of the algorithm is as follows [12]: System is initialised and we have t = t 0 , starting time, and x = x 0 , starting population, known.
2. Using a j as the propensity functions of the system evaluate a j (x) for each possible action j.
5. Replace t ← t + τ and x ← x + v j , where v j is the change of population when action j occurs.
6. Return to step 2 or end simulation.
In the case of the stochastic SIR model, where x = (S, I, R), we just replace a j with: a 0 (x) = β SI N a 1 (x) = γI (5) and v j with:

Markov Chain Monte Carlo
When trying to model an infectious disease with a stochastic SIR model the parameters of the model are going to be unknown. One of the ways we can estimate those parameters is Markov Chain Monte Carlo or MCMC -a class of computational methods -which in Bayesian Statistics are used to sample from the posterior distribution of some random variable. The term MCMC can be split into two parts, Monte Carlo and Markov Chain. Monte Carlo is another class of algorithms which are used to sample large amounts of random variables and get estimates of statistical properties for these random variables that resemble the real ones. Markov Chain is a stochastic process defined by a sequence or "chain" of events that happen in a set of time, which can be discrete or continuous. Markov Chains are characterised by the fact that they are "memory-less", meaning that the probability of being at a state, X n+1 , only depends on the previous state, X n , and not on the ones that came before it. Hence a discrete Markov Chain is defined by the following probability equation, also known as Markov "memory-less" property: Having defined the parts that comprise MCMC algorithms we can now describe how a typical MCMC algorithm works. MCMC algorithms are used to conduct bayesian inference so, given some realisation of a random variable X, they try to sample from the posterior distribution of a random variable Y . The Bayes theorem states that: where P (Y = y|X = x) is the posterior distribution of Y , P (X = x|Y = y) is the likelihood that data X are true for the given Y and P (Y = y) is the prior distribution of Y . If we define function p(y; x) = L(y; x)π(y), where L(y; x) = P (X = x|Y = y) and π(y) = P (Y = y), then we have up to a multiplication constant that does not depend on y, also known as the normalising constant.
It is very common that in many applications we may not know how to calculate the marginal probability distribution function of X, thus we cannot calculate the posterior distribution analytically. MCMC algorithms are used to construct intervals that are then used to estimate the posterior density via bayes rule and they have the advantage that you only need to know the target probability distribution function up to a normalising constant. MCMC algorithms work by generating samples θ t+1 from a proposal distribution h(·|θ t ) that depends only on the last sample generated (Markov "memory-less" property). Thus, they generate chains that after enough samples have been generated converge to the posterior distribution of θ (Monte Carlo simulation). In the specific MCMC algorithm we are going to be using, Metropolis-Hastings algorithm, a generated sample is accepted to the chain under a given probability, the acceptance probability, that depends on the proposal distribution h(·) and the target distribution P (·) which is proportional to the likelihood function times the prior for some given data X. If the generated sample of θ is not accepted, then the previously accepted sample is added to the posterior distribution instead.
Metropolis-Hastings algorithm [15] is one of the most common MCMC algorithms and is a generalisation of the Metropolis algorithm which is actually the first MCMC algorithm [19].
Algorithm 2 Metropolis-Hastings Algorithm 1. System is initialised and we have i = 0, starting time, and sample from prior distribution θ 0 ∼ v.
2. Using the proposal distribution we sample θ ∼ h(·|θ i ) and calculate the probability of 3. With probability p i we accept θ and θ i+1 ← θ and with probability 1 − p i we discard θ and θ i+1 ← θ i . 4. Increment i by 1 and return to step 2, or finish loop if a set number of samples from θ has been generated.
We can rewrite the probability of acceptance from the Metropolis-Hasting algorithm using the fact that the target distribution is proportional to the likelihood times the prior.
Thus, by Equation 9 we get that a simulated θ is more likely to be accepted the higher the posterior probability of θ given some data X compared to the posterior probability of the previously sampled θ i . Additionally, if the proposal distribution is symmetric, then h(θ i |θ ) = h(θ i |θ ) and the two terms cancel out. If the prior distribution is proportional to a constant, we also have v(θ ) = v(θ i ) and these two terms cancel out as well. Then the probability of acceptance depends only on the likelihood.
Though MCMC is a very useful set of algorithms it is hard to determine the proposal distribution and the prior distribution (initial conditions) to ensure convergence is achieved in a relatively small amount of steps and it is thus a computationally intensive process.

Approximate Bayesian Computation
Another way to sample from posterior distributions of parameters is Approximate Bayesian Computation or ABC [22]. ABC is a class of Bayesian computational methods used in parametric inference which, given a prior distribution of a set of parameters and some data, can estimate the posterior distribution with the advantage, compared to MCMC, that you do not need to know the posterior density up to a multiplication constant.
There are many algorithms included in the ABC class of algorithms, but the simplest one is the rejection sampling algorithm [4]. In the rejection sampling algorithm given some model g(·|θ) that simulates data depending on a parameter θ and a prior distribution for θ we draw n samples θ from the prior distribution and then use them to simulate data X i from g(·|θ i ) ∀ i ∈ [1, n]. If some selected summary statistics of X i are close enough to the summary statistics of the real data D, given a distance function, we accept θ i in the posterior distribution, otherwise we discard it. Thus, the algorithm requires summary statistics to compare on, a distance function to calculate how similar are the simulated and real data, as well as a threshold such that if distance between real and simulated data is less than threshold the sampled theta is accepted.
2. Sample X i ∼ g(·|θ i ) and then calculate S(X i ) for some summary statistics function S.
3. Calculate the distance between the summary statistics of the simulated and real data h(S(X i ), S(D)), where h is the distance function and D the observed data.
4. If the distance is less or equal than threshold , accept θ i otherwise do nothing. The estimated posterior distribution of θ from the rejection sampling algorithm will look like the following: where I is the indicator function. The exact posterior distribution of θ is θ ∼ v(·)f (S(D)|·). Thus, it is clear that when → 0 the approximated posterior distribution we get from ABC converges to the exact one.
In the case of the stochastic SIR model we can use a rejection sampling algorithm, by replacing g above with the Gillespie algorithm for SIR models, S with something like the number of infected and susceptible populations at each discrete time-point and h with the absolute difference of its arguments.

Sequential Monte Carlo
Sequential Monte Carlo are a set of Monte Carlo algorithms used to sample the posterior states of a Markov process [7]. Here we are going to use Sequential Monte Carlo, or Particle Filters as they are called, in the context of discrete time hidden Markov Chains. A hidden Markov Chain is described by a set of states X 0:n := (X 0 , ..., X n ), which are themselves a Markov Chain, but are not observed, and a set of states Y 0:n−1 := (Y 0 , ..., Y n−1 ), which are observed. Hidden Markov Models depend on the following distributions.
where µ is the initial distribution of X 0 , X p is sampled from a distribution that depends only on the previous state X p−1 and Y p is sampled from another distribution that depends only on the Markov Chain at time p which is X p . In the context of Hidden Markov Models some interesting statistics we want to retrieve are the likelihood of the observations Y p for p ∈ {0, n − 1} and get samples of the hidden Markov Chain X p for p ∈ {0, n}.
Concerning the likelihood of the observations there are ways to get the exact values analytically, but only for a few special cases. A more feasible way to estimate the likelihoods is using a Particle Filter ([6], [7]). A particle filter works by generating particles at each time-step from 0 to n that simulate X p , the hidden chain, and then calculate their weights according to how likely the particles are to match the actual X p given the observations Y p . Then in the next time-step only the particles with the highest weights are selected to sample from the next set of particles. The thinking behind particle filters is that at each time-step the system corrects itself by only choosing the best particles according to the observations. This can be less computationally intensive than other methods, since it requires less simulations to achieve the same or better results.
Algorithm 4 Particle Filter 1. Initiate the algorithm by setting Z 0 = 1, setting time t = 1 and sampling N particles 6. Return to step 2, or stop the algorithm if enough time-steps have passed.
Algorithm 4 describes a typical Particle Filter, where Z t = P (Y 0:t = y 0:t ) and we are interested in extracting Z n−1 , the likelihood of the observations under our model. We also see that Z n−1 is nothing more than the product of the mean of the weights at each time-step, where the weights are the probabilities of observing Y t if the particles describe the hidden Markov Chain X t . Then at each step, N indexes are sampled, which represent the index of the particles at the previous step, and the ones with the highest weights have the highest probability of being selected. The ones selected are used as the state of the hidden Markov Chain in the previous step and sampled from. Thus, only the particles that fit the data the best are selected and the rest are discarded.
Algorithm 5 Sampling a path 1. Set t = n − 1, sample a particle from the last time-step at random I n ∼ U (1, N ) and set X n ← ζ In n .

Follow back the path and set
, where A is the ancestral matrix, meaning A j t is the index of the ancestor particle of particle j at time t.

Set
4. Iterate t ← t − 1 and go to step 2, or end the algorithm if (X 0 , ..., X n ) has be sampled.
Then after Algorithm 4 has been ran, the Algorithm 5 can be run to extract a sample of the hidden Markov Chain X t . The algorithm is pretty straightforward, you select a random particle at the last time-step and then go backwards until you reach time 0 and have a completed chain.
There are various ways we can consider adapting our stochastic SIR model to a Hidden Markov Model. We can consider, for example, having X t be the hidden (S t , I t , R t ) values and Y t be the I t value, to simulate the scenario where only infections are observed. Another scenario could be observing only a fraction of the populations of each compartment and model that to retrieve the actual population numbers. Working with hidden Markov Chains makes the models more complicated and computationally expensive, but makes the epidemic models more realistic and, thus, more likely to be used in a real epidemic setting.

Particle Markov Chain Monte Carlo
A class of models used to estimate posterior distribution of random variables is Particle Markov Chain Monte Carlo or Particle-MCMC or PMCMC, which is a combination of Sequential Monte Carlo and an MCMC algorithm. It can be used to estimate the underlying parameters that influence the distributions of a hidden Markov Chain [7]. Therefore, a Hidden Markov Model's equations become the following, which is a variation of Equation 13.
where θ is unknown. A PMCMC algorithm estimates the parameter θ by using an MCMC algorithm to sample θ and using a Particle Filter to extract the likelihood of that θ. Hence, the Particle Filter is used to extract the likelihoods and a sample hidden Markov Chain at each time-step [2]. The rest of the algorithm is the same MCMC algorithm, as it would have been used on its own. Since we have a stochastic system which evolves through time, PMCMC algorithms take advantage of delaying prediction and using only the best trajectories and, thus, being less computationally intensive.
A common PMCMC algorithm is the particle marginal Metropolis-Hastings algorithm, which is the same as the Metropolis-Hastings algorithm with the addition of using the likelihoods extracted by a Particle Filter to accept a proposed θ . Define v as the prior distribution of parameter θ, then we have that the posterior distribution of θ is defined by: Consequently, the higher the posterior probability of θ , calculated with likelihood extracted from the particle filter, compared to the posterior probability of the previously accepted θ, the more likely to be accepted. 2. Using the proposal distribution we sample θ ∼ h(·|θ i ) and use θ to run the particle filter and extract Z N n,θ .
3. We calculate the probability of acceptance With probability p i we accept θ and θ i+1 ← θ and with probability 1 − p i we discard θ and θ i+1 ← θ i . 5. Increment i by 1 and return to step 2, until a set number of samples from θ has been generated.

Results
In this section I will provide examples of the uses of the algorithms mentioned in the "Methodology" Section 2 in an epidemiological setting. The models were fitted using simulated data, with known parameters, to validate their performance.
The code used to run the algorithms can be found here: https://github.com/GeorgeEfstathiadis/ Stochastic-Epidemic-Modelling

Noisy SIR model
I simulate a deterministic SIR epidemic with parameters β = 2 and γ = 1. The total population is 4820 with 4800 people in the susceptible compartment and 20 in the infectious compartment. The epidemic is simulated in continuous time, but only keeping the values at integer timepoints. I then proceed to add some Normal noise to imitate real world settings, where the disease is not following the exact deterministic SIR model. The noise is independent between different compartments and time-points and is described by e p ∼ N 3 (0, n ratio * X p * I 3 ), where p represents the time, X p the simulated, hidden, data at time p and n ratio is a predetermined ratio to enhance the noise depending on the size of the compartment at a given time. For this case I selected 0.01 as our n ratio . Thus the distribution g y (·|X p ) of the Hidden Markov Chain is N 3 (X p , 0.01 * X p * I 3 ). I proceed to model this simulated dataset using 2 different algorithms, ABC and PMCMC.

ABC
I assume β and γ have independent uniform priors in the range of (0, 5). For function g to simulate X i the Gillespie algorithm is used and for summary statistic function S the infected and removed counts at each time-point are used. Finally, for distance function the mean absolute difference function is used and for threshold , after a trial process, I ended up with 150. I gathered 1,000 samples to form the posterior distribution and the results were as follows.   Overall, the use of a non informative prior increased the computation time and we can clearly observe that it influenced the posterior distribution more than we would have wanted. For β the posterior mean is 2.16 with a 95% highest posterior density interval of (1.626, 2.768), while for γ the posterior mean is 1.12 with a 95% highest posterior density interval of (0.729, 1.52). The posterior median value of β is 2.12 and of γ is 1.10.
The algorithm achieves a decent performance, but this can be improved by using a more informative prior and a more complicated algorithm. In this setting ABC will not be able to scale well with more complicated scenarios.

PMCMC
I also model the noisy SIR dataset using a PMCMC algorithm. For the prior and proposal distribution of β and γ I use a multivariate normal distribution [23]. In consequence, the distribution looks similar to this: where Σ Σ Σ is the estimated posterior covariance matrix of β and γ and h is a constant that adjusts the covariance of the proposal distribution so that we have an acceptance rate between 10% and 25%. I also add a constraint to make sure our parameters are strictly larger than 0, while if any of the proposed parameters are below 0 we reject that sample.
I run the algorithm twice for around 1,000 epochs, one for estimating Σ Σ Σ and one for estimating h, without saving the results to the posterior distribution, a variation to the adaptive approach Haario and others formulated [14]. Hence, when I run the final algorithm, I have already estimated Σ Σ Σ, h and the prior distribution has as its mean the final sample of β and γ from the previous test run.
I then fix Σ Σ Σ and h and run the algorithm. I run the algorithm on 3 different chains to ensure convergence for 25,000 steps each. Since the process can be quite computationally intensive, I make sure to use parallel processing to receive our results faster. I also run tests on the number of particles needed to have consistent likelihoods returned from the particle filter for the same parameters. I noticed that difference in variance between the likelihoods for 10 and 100 particles was quite significant, but between 100 and 1000 was not significant. Thus I decided to run PMCMC with 100 particles in the interest of reducing the computational time.
After receiving the posterior samples of β and γ, I apply a burn-in of 1,000 samples, meaning I removed the first 1,000 samples received, in order to ensure that the parameters have reached convergence, before adding them to the posterior distribution. I also applied a thinning of 40 samples, meaning I only kept 1 sample every 40 samples for the posterior distribution, to reduce the auto-correlation between consecutive samples. After that I receive the estimated posterior distribution presented below. I also perform the Gelman-Rubin test to ensure all 3 chains have converged, which comes up with 1.00030734 and 1.00002084 for β and γ respectively, which are indeed sufficiently close to 1. I also calculate the effective sample size to see how many samples are independent by accounting for auto-correlation [17], which ends up being 1267 and 861 for each parameter.      The acceptance rate, before applying thinning is around 22%, while after applying thinning is around 100%. We also have the estimated posterior covariance is equal to: Σ Σ Σ = 0.00269008 0.0007161 0.0007161 0.00075565 (17) The posterior mean of β is 1.973 with a 95% highest posterior density interval of (1.884, 2.079), while for γ the posterior mean is 0.986 with a 95% highest posterior density interval of (0.927, 1.035). The posterior median of β is 1.972 and of γ is 0.984.

Under-reported SIR model
I use the same simulated dataset from the previous Subsection 3.1, with β equal to 2 and γ equal to 1 and with a population of 4820 with 4800 people in the susceptible compartment and 20 in the infectious compartment. Instead of adding Normal noise I only use a subset of our data to imitate real world under-reported data. I generate our observed data by applying a Binomial distribution to the, real, hidden data with a probability of observation being p obs . Thus, the distribution g y (·|X p ) of the Hidden Markov Chain is 3 independent Bin(X p , p obs ) distributions, one for each compartment. I proceed by modelling two different scenarios that might be useful in a real world epidemic. In the first model I assume that we know the probability of observing an infection and on the second scenario I assume that we don't know it and try to estimate the posterior of the probability of observation as well. For both scenarios I am using the PMCMC model and our observed data are generated using p obs equal to 0.1.

PMCMC with known missing rate
Similarly to the scenario discussed previously on noisy data, I use a multivariate normal distribution for the prior and proposal distributions of β and γ. I run two chains with around 1,000 epochs, one for estimating Σ Σ Σ and one for estimating h, without saving the results and choose an h and Σ Σ Σ than gives us a desirable acceptance rate as done in Subsection 3.1.2.
After some analysis of the likelihoods I decide to use 1,000 particles and run the algorithm on 3 different chains with 25,000 steps each to ensure convergence. When the algorithm ends, I apply a burn-in of 1,000 samples and a thinning of 40. I also perform the Gelman-Rubin test diagnostic to ensure convergence on the 3 chains and I get 1.00030734 and 1.00002084 statistics for β and γ respectively, which are sufficiently close to 1. I also calculate the effective sample size, which ends up being 1050, 1300 for each parameter.      Before applying the thinning we have an acceptance rate of around 12% and after applying thinning we get an acceptance rate of around 95%. We also have the estimated posterior covariance being: Σ Σ Σ = 0.00653857 0.00305867 0.00305867 0.0019865 (18) For β a posterior mean of around 2.036 is achieved with a 95% highest posterior density interval of (1.907, 2.205), while for γ a posterior mean of 1.008 is achieved with a 95% highest posterior density interval of (0.93, 1.094). The posterior median of β is around 2.036 and the posterior median of γ is around 1.003.

PMCMC with unknown missing rate
Modelling the scenario with an unknown missing rate has one difference from the previous scenario; there is one extra variable's posterior distribution to estimate. I follow the same steps as in the previous subsections but the prior and proposal distributions are now 3-dimensional since they include p obs as the third variable.
After some analysis of the likelihoods I decide to use 100 particles again and run the algorithm on 3 different chains with 25,000 steps each to ensure convergence. When the algorithm ends, I apply a burn-in of 1,000 samples and a thinning of 40. I also perform the Gelman-Rubin test to ensure convergence and end up with 1.00040369, 1.00070195 and 1.00025463 for β, γ and p obs respectively, which are all sufficiently close to 1. I also calculate the effective sample size, which ends up being 1275, 1732 and 2622 for each parameter.        Before thinning the acceptance rate is around 23%, while after applying thinning the acceptance rate is around 86%. Meanwhile, the estimated posterior covariance is equal to: The posterior mean of β is 2.034 with a 95% highest posterior density interval of (1.79, 2.311), the posterior mean of γ is 1.006 with a 95% highest posterior density interval of (0.877, 1.16) and the posterior mean of p obs is 0.101 with a 95% highest posterior density interval of (0.098, 0.105). The posterior median of β is 2.018, of γ is 0.977 and of p obs is 0.101.

Under-reported SEIR model
I am now going to use a slightly more complicated model than the SIR, which was discussed in Subsection 2.1, the SEIR model. This model includes an extra parameter, α, which signifies the rate at which people from the susceptible group enter the exposed group. I will work with under-reported data again, generated by applying the binomial distribution to the simulated data, but I am only going to model the scenario where the probability of observation is known. This is done as proof of concept and expanding it to the scenario where p obs is unknown would be modelled the same way.
I generate our data from a deterministic SEIR model with β = 4, γ = 1 and α = 1. Total population is 4820 people, with 4800 people in the susceptible group, 20 in the infected group and 0 in both the exposed and removed group. I then apply a binomial distribution with p obs = 0.1 to imitate real world conditions, when data are under-reported. Thus, the distribution of the hidden Markov Chain is the same as in the previous section.
I model this scenario using the PMCMC algorithm. For the prior and proposal distribution of β, γ and α I use a multivariate normal distribution, as in the previous subsections, which is 3-dimensional. I run the algorithm with 3 chains using 100 particles and for 25,000 steps each. After the algorithm is finished I apply a burn-in of 1,000 samples and thinning of 20. I also perform the Gelman-Rubin test to ensure convergence, which ends up being 1.00103306, 1.00010529 and 1.00098897 for β, γ and α respectively, which are sufficiently close to 1. I also calculate the effective sample size, which ends up being 752, 2065 and 835 for each parameter.
For β a posterior mean of 4.294 is achieved with a 95% highest posterior density interval of (3.462, 5.443), for γ a posterior mean of 0.946 is achieved with a 95% highest posterior density interval of (0.803, 1.091) and for α a posterior mean of 1.064 is achieved with a 95% highest posterior density interval of (0.797, 1.362). The posterior medians are 4.178, 0.948 and 1.034 for β, γ and α respectively.

Experiments on information reduction
The following section examines various experiments that were performed on reducing the amount of data the algorithms are trained on to pinpoint the effects of these changes on convergence. It is of interest to find out how much the observed datasets can be modified, with the chains still converging to their posterior densities.
In order to assess the performance of the algorithm in the modified datasets I used the density plots of the posterior density and the trajectory plots. I also calculated the posterior means and high posterior density intervals of the underlying parameters, as well as their posterior mean squared error defined by: where θ i is the i th true parameter used and θ i,j it the j th posterior sample of the i th parameter from the estimated posterior samples.
For the purposes of this experimentation on the amount of information needed I am using the adaptive version of the algorithm used in the previous subsections proposed by Haario and others [14]. In essence, I fix h at a sensible value, using prior knowledge, and I adapt the posterior covariance of the parameters Σ, while the algorithm runs. Thus, the formula for Σ now is: where I n is the identity matrix with n dimensions, n being the dimensions of θ. is fixed to be 0.0001 and t 0 to be 1,000. Thus, we do not start adapting the posterior covariance until the first 1,000 steps have passed and then applying a burn-in of 1,000 steps to each chain.
Hence, the new proposal distribution inside the PMCMC will be θ t+1 ∼ N n (θ t , hΣ Σ Σ t ). The difference with the results we got from the previous subsections is that now the results are no longer Markov Chains, since a sample θ t+1 will depend on all previous samples, due to the adaptive nature of the posterior covariance.

Noise level
The results of increasing the noise, for the same SIR noisy model discussed above, are investigated here. Originally, convergence was achieved on the dataset with a noise level of 0.01. For this experiment, I ran the algorithm on datasets created using a n ratio of 0.05, 0.1, 0.15, 0.2 and 0.25. I run one chain for each n ratio used, with 6,000 samples each with a burn-in of 1,000 samples and a thinning of 10 samples for the visualisations only.
For n ratio equal to 0.05, the posterior mean of β was 2.018 with a high posterior density interval of (1.917, 2.142), while the posterior mean of γ was 1.009 with a high posterior density interval of (0.954, 1.068). The posterior mean squared error of β and γ were 0.004008 and 0.00097 respectively.         For all the different noise levels the chains were able to produce fairly good results. While the posterior means remain close to the true value and the high posterior density intervals still include the true values, we can observe that as the noise level increases, so does the posterior variance of the parameters. We can see that visually from the density plots and the credible interval of the trajectories, which get wider. We can also deduce that from the posterior mean squared error increasing. The posterior bias squared remains relatively stable, thus the posterior variance is causing the increase. This is to be expected, since as the noise increases the more unsure we are about the posterior estimates. Thus, while increasing the noise, the posterior estimates are still very close to the true values, the posterior mean squared error remains relatively low and the sampled trajectories are still close to the truth.

Probability of observation level
In the following subsection the results on experimenting with the probability of observation, p obs , for the under-reported SIR model with known p obs are presented. In the previous subsection the algorithm was fitted for a p obs equal to 0.1. Now, I will investigate what happens if p obs is equal to 0.075, 0.05, 0.025 and 0.01. One chain is run for each p obs used, with 6,000 samples each with a burn-in of 1,000 samples and a thinning of 10 samples for the visualisations only.
For the 0.075 p obs , the posterior mean of β was 2.013 with a high posterior density interval of (1.84, 2.217), while the posterior mean of γ was 1.0 with a high posterior density interval of (0.883, 1.104). The posterior mean squared error of β and γ were 0.009624 and 0.003027 respectively.        From the results above we can deduce that while decreasing the probability of observation in the under-reported dataset, we still achieve good results. Similarly to the tests run in Subsection 3.4.1 on increasing the noise, the posterior variance of the results seems to increase, while the posterior means remain very close to their true values. The posterior mean squared errors seem to generally be increasing when we reduce p obs and the variability of the sampled parameters and trajectories increases, as it can be seen from their credible intervals. This can be attributed to the reduction of information and data-points read by the model every time we reduce the probability of observation. Yet, fairly good results are still achieved, even with this limited amount of data.

Time-points observed level
The effect of reducing the time-points observed in our predictions is additionally investigated in this subsection. In epidemiology it is essential to gather information about the disease as early as possible, so using a simple SIR model with noisy data I try to find at which point during the disease we would be able to make accurate predictions about the disease's dynamics. I use the same dataset as described in Subsection 3.1, which had n ratio of 0.01, and I test reducing the time-points from 15 originally to 11, 7 and 3. I run one chain for each case with 6,000 steps each with a burn-in of 1,000 samples and a thinning of 10 samples for the visualisations only.
When using the 11 first time-points, the posterior mean of β was 1.948 with a high posterior density interval of (1.441, 2.504), while the posterior mean of γ was 0.958 with a high posterior density interval of (0.521, 1.424). The posterior mean squared error of β and γ were 0.082036 and 0.063657 respectively.      We deduce that when using the 11 first time-points we still get decent results, although the posterior variance increases significantly and the credible intervals of the sampled trajectories get much wider, and that could be since the peak of the disease has passed by this point and it is still possible to infer the parameters' posterior distributions from the data available. When reducing the time-points to 7, we notice that there is another increase in the posterior covariance of the parameters and the parameters are being overestimated. This could be because of the very limited amount of time-points. When only using the first 3 time-points, the posterior estimates of the parameters are overestimated and are very far from the true values, while their posterior variances are significantly increased, to a point that it seems as if the model cannot possibly converge to the truth. This is probably due to the lack of identifiability in the dataset -for this amount of information there are possibly many different parameters that could explain the observed data. That could be due the limited time-points and specifically for not having the peak of the disease included in the data provided, meaning that in general the data could result into very different epidemics.

Discussion
In this project the use of Particle Markov Chain Monte Carlo in stochastic epidemic modelling is investigated, using multiple simulated datasets and hidden Markov Chains for different possible scenarios. Additionally, an adaptive approach to modelling the posterior covariance of the parameters is explored in order to model scenarios where the information in the data read by the model would sequentially reduce. The algorithm was able to retrieve the estimated posterior densities, even with a very limited amount of information on the data. For models with few parameters very good posterior densities were received, centered at the true values, and credible intervals around the trajectories that were very accurate. When reducing the amount of data observed through the hidden Markov Chains, the uncertainty of the results was evidently increasing. Yet, very accurate posterior estimates were still received, when the identifiability requirement was not violated. These results are promising, since in a real world setting the uncertainty of the observed cases of a disease is expected to be very close to the ones modelled here. Thus, any parties interested in getting some of the characteristics underlying an epidemic's dynamics could use this algorithm to reach some estimates quickly, with a minimum amount of data that are accessible to anyone.
While the algorithm can be very useful in certain settings, it is also not without limitations. One such limitation is that the algorithm is very computationally intensive. MCMC by nature is a class of algorithms that due to Monte Carlo, requires generating a large amount of samples and when extending it to include the particle filter, which generates a large amount of particles at each step, it requires significant computational resources. In this project only small datasets were simulated and simple models used, but it was still difficult to get an acceptance rate of more than 20%, due to the complicated nature of the epidemic models. Thus, it requires an even larger sample than those used in simple MCMC algorithms to achieve convergence and achieve an uncorrelated estimate of the posterior density. In addition, increasing the dimensions of the parameters makes the algorithm even more computationally intensive. During this project it was also attempted to fit a model with multiple subgroups, which had different rates of contact and infection between the subgroups. It proved to be challenging to deliver good estimates of the posterior covariance and achieve convergence, so due to time constraints the experiment was dropped. This could be due to the use of the Metropolis-Hastings algorithm and perhaps a different MCMC algorithm would be better suited for a large number of parameters. In recent years new MCMC algorithms have been proposed that are better suited for multi-dimensional inference [9] and for ODE based models [24] which could be an interesting extension to the work presented here by implementing the particle filter in combination with those algorithms.
In relation to possible extensions to this work, a couple of topics were aimed to be included in this work, but this has not been feasible due to the limited time-frame. It would be interesting to try to fit an SIR like model with the model trained only in a single compartment. This should probably be the infectious compartment, which could be very useful, since in real world settings the infectious count data are sometimes the only available. In addition, it would be intriguing to try to fit a model with multiple subgroups, but with observing only the totals by compartment, to reduce further the amount of information needed. Possible extensions that could be of interest, would be to try to fit more complicated models, like the SIS [20] or SEIRS [3], or models with changing populations due to demographics [18]. These scenarios would add a layer of complexion due to the infectious curve having multiple modes and being harder to model, which is necessary when dealing with diseases that have an element of seasonality, like the influenza.