Markov chains as a simulation technique for epidemic growth Karolina K˛esik Institute of Mathematics Silesian University of Technology Kaszubska 23, 44-100 Gliwice, Poland Email: karola.ksk@gmail.com Abstract—Modeling various phenomena occurring in nature [9]–[11], whether the deployment of service stations in some allows us to predict future effects in reality. One of such area due to existing roads and their traffic [12]. example is modeling the growth of infection in a given Another application is artificial neural networks, where population depending on various parameters. In this article, we show the use of discrete Markov chains in order to model the there is a specific variation thereof with a stochastic note. epidemic with distinction into four states in which individuals In [13], [14], the authors model such networks and analyze in the population may be. More accurately – healthy, infected, the flow of information. Hybrids are also created, such as the sick and recovered. The article presents a mathematical model connection of heuristics with neural networks [15]. Not only describing the phenomenon together with a calculation example neural networks, but other artificial intelligence algorithms are and simulations. The obtained results were described and discussed. used in practical terms, an example of which is optimization green computing awareness [16]. Another interesting idea Keywords: discrete Markov chains, epidemic growth, mathematical is modeling of contaminant transport in porous media and model, simulation monitoring of water quality [17], [18]. In this paper, we describe one of the classic stochastic I. I NTRODUCTION tools such as Markov chains. We show their use in modeling Stochastic processes are a family of random variables de- the epidemic phenomenon and pay attention to their use in fined in a certain probabilistic space. It is the field of the forecasting future phenomena. probability theory, which in today’s science has become one of II. M ARKOV CHAINS the most dynamically developed fields due to the numerous ap- plications in optimization techniques or artificial intelligence. The process {Xn , n ≥ 0} of state space S is called Such processes enable defining various phenomena based on a the discrete Markov chain, if for each n ∈ {0, 1, . . . }, the certain probability of its occurrence, as well as its components. following equation occurs The field of stochastics found its place in the theory of P {Xn+1 = in+1 |Xn = in , . . . , X0 = i0 } optimization, where Lévy and Poisson processes were used (1) = P {Xn+1 = in+1 |Xn = in } in modeling various types of phenomena occurring in nature. Optimization theory has gained the most, where techniques for all possible states i0 , . . . , in+1 ∈ S. inspired by nature have been created to this day. An example It is a mathematical model of a random phenomenon is the cuckoo algorithm intended originally to search for ex- evolving over time in such a way that the past affects the tremes of function. The movement of these birds was modeled future only through the present. This model has state space S, using the Lévy flights. In [1]–[3], the author presented the where we can give the following properties idea of searching key-points on images which can be used 1) S it is a finite or at most a countable set of states, for feature extractions. Proposed idea is important in security 2) S = {0, 1, 2, 3, . . .}, because using stochastic algorithms almost guarantees finding 3) Xn = i, what means, that at time n, the process is in other features due to the random placement of individuals the state of i. in the population at the beginning of the algorithm. With Let us define the initial distribution of states, i.e. the dis- similar motives in [4], [5], the author showed the use of tribution of the variable X0 . A probability vector as π = these algorithms in two-dimensional games, where such an [p0 , p1 , . . . ]. For each state i ∈ S, we have algorithm can be used to model the opponent’s movement. Additionally, in [6]–[8], heuristic approach was used in dif- pi = P {X0 = i}. (2) ferent games. Similar algorithms are genetic ones based on chromosomes. Their use is visible in route planning programs In addition, we define the probability of transition between two specific states. It is determined using a matrix P = (pij ). c 2019 for this paper by its authors. Use permitted under Creative If S is a set defined as (1, 2, . . . , m), than a matrix P has Commons License Attribution 4.0 International (CC BY 4.0) dimension equal m × m. As pij , we can define the probability 1 of going from one state i to another j in one step. Generalizing IV. E XPERIMENTS this, we have A. Numerical example pij = P {Xn+1 = j|Xn = i}. (3) Let us assume that at the beginning of our experiment, the Using the initial vector π and the probability matrix P, we population is made up of a thousand people, where 495 people define the Markov chain. And its behavior is defined using a are healthy, 3 are infected and 2 are sick. Hence, the initial stochastic matrix P = (pij ). vector is  We assume that every element of such a matrix must be P0 = 495 3 2 0 . (8) greater than 0 and the sum of elements of each row must be The stochastic matrix is defined as follows 1. This matrix has the following form     0.8 0.1 0.1 0 p11 p12 . . . p1m  0 0.35 0.4 0.25   p21 p22 . . . p2m  P=  0 . (9) 0 0.85 0.15  P= . (4)   .. .. ..   .. . . .  0 0.05 0.02 0.93 pm1 pm2 . . . pmm After the first iteration, there is If the Markov chain is in the state i at n, then the probability  0.8 0.1 0.1 0  of being in the state j after k periods is defined as follows  0 0.35 0.4 0.25  P0 P = 495 3 2 0   0 = P (Xn+k = j|Xn = i) = P (Xk = j|X0 = i) = pij (k), (5) 0 0.85 0.15  0 0.05 0.02 0.93 where probabilities are independent of n, pij (k) is an element  at position (i, j) in Pk . = 396 51 52 1 . (10) III. M ARKOV CHAINS FOR EPIDEMIC SIMULATIONS After the second iteration, there is The spread of infection can be represented using a math-   ematical model – a discrete Markov chain. In this model, 0.8 0.1 0.1 0 we divide the population into four groups marked as states   0 0.35 0.4 0.25  396 51 52 1   0 = S={0, 1, 2, 3}, i.e. 0 0.85 0.15  • 0 – susceptible sus, 0 0.05 0.02 0.93  • 1 – infected inf , = 317 57 104 21 . • 2 – sick sick, (11) • 3 – recovered rec. B. Simulation experiments Unfortunately, there is no way to create a perfect model, that is why we create several assumptions that simplify this Simulations were conducted for two stochastic matrices. model The first of these was described by the Eq. (9) and the second • a healthy person can become infected with an infection, one has the following form   • an individual in the infected group may go only to the 0.8 0.1 0.1 0 disease state,  0 0.59 0.4 0.01  P= . (12) • a sick subject can leave a group of patients only through  0 0 0.99 0.01  complete recovery, 0 0.05 0.02 0.93 • recovery guarantees immunity, Both matrices are similar but differ in the selected values. • immunity is not inherited In the second matrix, the probability of transition from the • age, sex and social status do not affect the probability of infected to recovered and sick to the recovered state drastically infection, changed, which means minimal chance of recovery. In all • climate or demographic changes do not affect the epi- simulations, we assumed that the population is composed of demic. 500 individuals. The effect of the difference in people needed For such assumptions, the stochastic matrix can be defined as   for the spread of the disease on the entire population was psus,sus psus,inf psus,sick psus,rec examined. Tests were performed for the population composed  pinf,sus pinf,inf pinf,sick pinf,rec  of 1000000 individuals in population and 1 sick for two   psick,sus psick,inf psic,sick psick,rec   (6) different step (time) values – {25, 150}. The results are shown prec,sus prec,inf prec,sick prec,rec in the form of diagrams in the Fig. 1–4. It is easy to see where p is the probability of changing states at t to t + 1. An that the charts are identical regardless of the initial vector initial vector can be represented as parameters (for the same probability matrix). Hence, the  simple conclusion that Markov models allow simulation of the P0 = nsus ninf nsick nrec , (7) epidemic phenomenon, but the initial vector has little effect. where n is the number of individuals in specific group at the Mainly stochastic matrices play a role, thus estimating the initial stage. probability of transition between selected states. 2 Figure 1: Measurements for the first matrix and population {495, 3, 2, 0}. Figure 2: Measurements for the second matrix and population {495, 3, 2, 0}. Figure 3: Measurements for the first matrix and population {1000000, 0, 1, 0}. Figure 4: Measurements for the second matrix and population {1000000, 0, 1, 0}. Note that in the probability of transition from infected and infection can evolve, and then these matrices can change. sick to healed state are 25% and 15%. In the population the Unfortunately, the classic approach to Markov’s chains does average number of healthy is around 80% people and it stays not modify the matrix during operation, although there are independent of the step. Large differences in jumps between models that can do it. Impact on these matrices will result in the population can be seen in the first 20 steps, where the a much better realignment of the model. population is infected and gets ill. The further step, the more For the experiments we have carried out, we have performed stable the chart is, which may be due to the lack of changes statistical tests. On the significance level α = 0.1, we also in the stochastic matrix. From a practical point of view, the verify the hypothesis about the compatibility of the data distri- 3 Table I: Statistical tests for susceptible table states can be huge, but then the problem arises to create such Statistic p-value a matrix and find its coefficients. The simulations showed that Anderson-darling 12173.87 0 there is a large impact of even a small change in the main Cramer-von Mises 34.27 0 Kolmogorov-Smirnov 0.84 1.59 · 10−95 matrix. Kuiper 0.92 6.11 · 10−114 Pearson χ2 1573.91 5.65 · 10−328 R EFERENCES Watson U 2 11.52 0 [1] D. Połap and M. Woźniak, “Voice recognition by neuro-heuristic method,” Tsinghua Science and Technology, vol. 24, no. 1, pp. 9–17, Table II: Statistical tests for infected table 2019. [2] D. Połap, “Model of identity verification support system based on voice Statistic p-value and image samples,” j jucs, vol. 24, no. 4, pp. 460–474, 2018. Anderson-darling 406.59 0 [3] S. Coco, A. Laudani, F. Riganti Fulginei, and A. Salvini, “Team Cramer-von Mises 45.77 1.22 · 10−15 problem 22 approached by a hybrid artificial life method,” COMPEL- Kolmogorov-Smirnov 0.93 8.45 · 10−117 The international journal for computation and mathematics in electrical Kuiper 0.93 3.19 · 10−115 and electronic engineering, vol. 31, no. 3, pp. 816–826, 2012. Pearson χ2 2084.21 4.73 · 10−438 [4] A. Winnicka, “The opponent’s movement mechanism in simple games Watson U 2 11.69 0 using heuristic method,” Symposium for Young Scientists in Technology, Engineering and Mathematics (SYSTEM 2018). [5] D. Połap, M. Wozniak, C. Napoli, and E. Tramontana, “Is swarm intelligence able to create mazes?” International Journal of Electronics bution with the Gamma distribution with the shape parameter and Telecommunications, vol. 61, no. 4, pp. 305–310, 2015. equal to 1 and the scale parameter equal to 2, which results are [6] V. V. Zakharov and V. A. Shirokikh, “Heuristic evaluation of the char- acteristic function in the cooperative inventory routing game,” Journal presented in Tab. I–IV. According to the tests carried out, in the on Vehicle Routing Algorithms, vol. 1, no. 1, pp. 19–32, 2018. case of each table presenting the process of healthy, infected, [7] A. Iordan, “A comparative study of the a* heuristic search algorithm sick, recovered, at the significance level α = 0.1, we can reject used to solve efficiently a puzzle game,” in IOP Conference Series: Materials Science and Engineering, vol. 294, no. 1. IOP Publishing, the hypothesis about the compatibility of the distribution of 2018, p. 012049. data with the Gamma distribution with the parameter equal to [8] D. Połap, M. Woźniak, C. Napoli, E. Tramontana, and R. Damaševičius, 1 and scale parameter 2. “Is the colony of ants able to recognize graphic objects?” in International Conference on Information and Software Technologies. Springer, 2015, pp. 376–387. V. C ONCLUSION [9] V. Roberge, M. Tarbouchi, and G. Labonté, “Fast genetic algorithm path planner for fixed-wing military uav using gpu,” IEEE Transactions on Discrete Markov chains allow us to model phenomena Aerospace and Electronic Systems, 2018. occurring in nature. Each step depends on the previous one, [10] S. Yang, H. Wang, Z. Ren, S. Mu, J. Li, and J. Zhao, “Distribution although there is no change in probabilities during operation. network inspection route planning and the application of inspection automatic control technique,” in 2018 IEEE 3rd Advanced Information This is a quite serious shortcoming in predicting the future Technology, Electronic and Automation Control Conference (IAEAC). effects of the model. In the analyzed case of epidemic spread, IEEE, 2018, pp. 1312–1315. such action resulted in the absence of a possible mutation of [11] C. Napoli and E. Tramontana, “A cloud oriented support for motion de- tection in video surveillance systems using computational intelligence,” the disease. However, it is worth noting that if the disease in International Conference on Information and Software Technologies. started to be fatal after a certain amount of time, in the case Springer, 2016, pp. 453–463. of the second matrix, it could kill almost the entire population, [12] T. H. Tran, G. Nagy, T. B. T. Nguyen, and N. A. Wassan, “An efficient heuristic algorithm for the alternative-fuel station location problem,” as opposed to the first one. European Journal of Operational Research, vol. 269, no. 1, pp. 159–170, This type of modeling of phenomena may allow us to 2018. improve the predictions of some phenomena that depend only [13] P. Baskar, S. Padmanabhan, and M. Syed Ali, “Novel delay-dependent stability condition for mixed delayed stochastic neural networks with on selected states and the table of probabilities. The number of leakage delay signals,” International Journal of Computer Mathematics, pp. 1–14, 2018. [14] Z. Liu, “Design of nonlinear optimal control for chaotic synchronization Table III: Statistical tests for sick table of coupled stochastic neural networks via hamilton–jacobi–bellman equation,” Neural Networks, vol. 99, pp. 166–177, 2018. Statistic p-value [15] S. Chatterjee, S. Sarkar, N. Dey, A. S. Ashour, and S. Sen, “Hybrid Anderson-darling 22104.23 0 non-dominated sorting genetic algorithm: Ii-neural network approach,” Cramer-von Mises 49.74 0 in Advancements in Applied Metaheuristic Computing. IGI Global, Kolmogorov-Smirnov 0.99 3.04 · 10−132 2018, pp. 264–286. Kuiper 0.99 9.11 · 10−131 [16] E. Okewu, S. Misra, R. Maskeliūnas, R. Damaševičius, and Pearson χ2 147.19 2.61 · 10−24 L. Fernandez-Sanz, “Optimizing green computing awareness for envi- Watson U 2 12.35 0 ronmental sustainability and economic security as a stochastic optimiza- tion problem,” Sustainability, vol. 9, no. 10, p. 1857, 2017. Table IV: Statistical tests for recovered table [17] V. Nourani, S. Mousavi, D. Dabrowska, and F. Sadikoglu, “Conjunction of radial basis function interpolator and artificial intelligence models for Statistic p-value time-space modeling of contaminant transport in porous media,” Journal Anderson-darling 3601.29 0 of hydrology, vol. 548, pp. 569–587, 2017. Cramer-von Mises 46.94 6.66 · 10−16 [18] D. Dabrowska, ˛ A. Witkowski, and M. Sołtysiak, “Representativeness of Kolmogorov-Smirnov 0.95 3.19 · 10−122 the groundwater monitoring results in the context of its methodology: Kuiper 0.95 2.19 · 10−122 case study of a municipal landfill complex in poland,” Environmental Pearson χ2 1967.18 8.61 · 10−413 Earth Sciences, vol. 77, no. 7, p. 266, 2018. Watson U 2 12.04 0 4