Composite Monte Carlo decision making under high uncertainty of novel coronavirus epidemic using hybridized deep learning and fuzzy rule induction.

In the advent of the novel coronavirus epidemic since December 2019, governments and authorities have been struggling to make critical decisions under high uncertainty at their best efforts. In computer science, this represents a typical problem of machine learning over incomplete or limited data in early epidemic Composite Monte-Carlo (CMC) simulation is a forecasting method which extrapolates available data which are broken down from multiple correlated/casual micro-data sources into many possible future outcomes by drawing random samples from some probability distributions. For instance, the overall trend and propagation of the infested cases in China are influenced by the temporal-spatial data of the nearby cities around the Wuhan city (where the virus is originated from), in terms of the population density, travel mobility, medical resources such as hospital beds and the timeliness of quarantine control in each city etc. Hence a CMC is reliable only up to the closeness of the underlying statistical distribution of a CMC, that is supposed to represent the behaviour of the future events, and the correctness of the composite data relationships. In this paper, a case study of using CMC that is enhanced by deep learning network and fuzzy rule induction for gaining better stochastic insights about the epidemic development is experimented. Instead of applying simplistic and uniform assumptions for a MC which is a common practice, a deep learning-based CMC is used in conjunction of fuzzy rule induction techniques. As a result, decision makers are benefited from a better fitted MC outputs complemented by min-max rules that foretell about the extreme ranges of future possibilities with respect to the epidemic.

Introduction

On top of devastating health effects, epidemic impacted hugely on world economy. In the Ebola outbreak between 2014–2016 where more than 28,000 and cases were suspected and 10,000 deaths in West Africa [1], $2.2 billion was lost [2]. On the other hand, SARS took over 648 lives from China including Hong Kong and 700 lives worldwide between 2002 and 2003 [3]. Its losses on global economy are up to a huge $100 billion, 1% and 0.5% dips of GDPs in Chinese and Asian domestic markets respectively [3]. Although the current coronavirus (codename: NCP or COVID-19) epidemic is not over yet, its economy impact is anticipated by economists from IHS Market to be far worse than that of SARS outbreak in 2003 [4]. The impact is so profound that will lead to factories shut down, enterprises bankruptcy especially those in tourism, retail and F&B industries, and suspensions or withdrawals in long-term investment, if the outbreak cannot be contained in time. Since the first case in December 2019, the suspected cases and deaths around the world skyrocketed to over 76395 confirmed cases and 2348 deaths, mostly in China, by the time of writing this article.

An early intervention measure in public health to thwart the outbreak of COVID-19 is absolutely imperative. According to a latest mathematical model that was reported in research article by The Lancet [5], the growth of the epidemic spreading rate will ease down if the transmission rate of the new contagious disease can be lowered by 0.25. Knowing that the early ending the virus epidemic or even the reduction in the transmission rate between human to human, all governments especially China where Wuhan is the epicentre are taking up all the necessary preventive measures and all the national efforts to halt the spread. How much input is really necessary? It is fundamentally challenging in making a sound decision. When an epidemic has just started, available data is scarce, information is limited, and uncertainty is high. Unlike common Influenza, there is neither precedent case exactly the same nor sufficient knowledge. Many decision makers could only take references from SARS which is by far the most similar virus to COVID-19. However, it is difficult as the characteristics of the virus are not fully known, its details and about how it spreads are gradually unfolded from day to day. Given the limited information on hand about the new virus, and the ever evolving of the epidemic situations both geographically and temporally, it boils down to grand data analytics challenge this analysis question: how much resources shall be enough to slow down the transmission? This is a composite problem that requires cooperation from multi-prong measures such as medical provision, suspension of schools, factories and office, minimizing human gathering, limiting travel, strict city surveillance and enforced quarantines and isolations in large scales. There is no easy single equation that could tell the amount of resources in terms of monetary values, manpower and other intangible usage of infrastructure; at the same time there exist too many uncertain variables from both societal factors and the new development of the virus itself. For example, the effective incubation period of the new virus was found to be longer than a week, only some time later after the outbreak. Time is an essence in stopping the epidemic so to reduce its damages as soon as possible. However, uncertainties are the largest obstacle to obtain an accurate model for forecasting the future behaviours of the epidemic should intervention apply. In general, there is a choice of using deterministic or stochastic modelling for data scientists; the former technique based solely on past events which are already known for sure, e.g. if we know the height and weight of a person, we know his body mass index. Should any updates on the two dependent variables, the BMI will be changed to a new value which remains the same for sure no matter how many times the calculation is repeated. The latter is called probabilistic or stochastic model — instead of generating a single and absolute result, a stochastic model outputs a collection of possible outcomes which may happen under some probabilities and conditions.

Background

Deterministic model is useful when the conditions of the experiment are assumed rigid. It is useful to obtain direct forecasting result from a relatively simple and stable situation in which its variables are unlikely to deviate in the future. Otherwise, for a non-deterministic model, which is sometimes referred as probabilistic or stochastic, the conditions of a future situation under which the experiment will be observed, are simulated to some probabilistic behaviour of the future observable outcome. For an example of epidemic, we want to determine how many lives could be saved from people who are infected by a new virus as a composite result of multi-prong efforts that are put into the medical resources, logistics, infrastructure, spread prevention, and others; at the same time, other contributing factors also matter, such as the percentage of high-risk patients who are residing in that particular city, the population and its mobility, as well as the severity and efficacy of the virus itself and its vaccine respectively. Real-time tools like CDC data reporting and national big data centres are available with which any latest case that occurs can be recorded. However, behind all these records are sequences of factors associated with high uncertainty. For example, the disease transmission rate depends on uncertain variables ranging from macro-scale of weather and economy of the city in a particular season, to the individual’s personal hygiene and the social interaction of commuters as a whole. They are dynamic in nature that change quickly from time to time, person to person, culture to culture and place to place. The phenomena can hardly converge to a deterministic model. Rather, a probabilistic model can capture more accurately the behaviours of the phenomena. So for epidemic forecast, a deterministic model such as trending is often used to the physical considerations to predict an almost accurate outcome, whereas in a non-deterministic model we use those considerations to predict more of a probable outcome that is probability distribution oriented.

In order to capture and model such dynamic epidemic recovery behaviours, stochastic methods ingest a collection of input variables that have complex dependencies on multiple risk factors. The epidemic recovery can be viewed in abstract as a bipolar force between the number of populations who has contracted the disease and the number of patients who are cured from the disease. Each group of the newly infested and eventually cured (or unfortunately deceased) individuals are depending on complex societal and physiological factors as well as preventive measures and contagious control. Each of these factors have their underlying and dependent factors carrying uncertain levels of risks. A popular probabilistic approach for modelling​ the complex conditions is known as Monte Carlo (MC) simulation which provides a means of estimating the outcome of complex functions by simulating multiple random paths of the underlying risk factors. Rather than deterministic analytic computation, MC uses random number generation to generate random samples of input trials to explore the behaviour of a complex epidemic situation for decision support. MC is particularly suitable for modelling epidemic especially new and unknown disease like COVID-19 because the data about the epidemic collected on hand in the early stage are bound to change. In MC, data distributions are entered as input, since precise values are either unknown or uncertain. Output of MC is also in a form of distribution specifying a range of possible values (or outcome) each of which has its own probability at which it may occur. Compared to deterministic approach where precise numbers are loaded as input and precise number is computed to be output, MC simulates a broad spectrum of possible outcomes for subsequent expert evaluation in a decision-making process.

Monte Carlo Simulation for epidemics

Recently as epidemic is drawing global concern and costing hugely on public health and world economy, the use of MC in epidemic modelling forecast has become popular. It offers decision makers an extra dimension of probability information so called risk factors for analysing the possibilities and their associated risk as a whole. Decades ago, there has been a growing research interest in using MC for quantitatively modelling epidemic behaviours. Since 1957, Bailey et al. was among the pioneers in formulating mathematical theory of epidemics. Subsequently in millennium, Andersson and Britton [6] adopted MC simulation techniques to study the behaviour of stochastic epidemic models, observing their statistical characteristics. In 2003, House et al. attempted to estimate how big the final size of an epidemic is likely to be, by using MC to simulate the course of a stochastic epidemic. As a result, the probability mass function of the final number of infections is calculated by drawing random samples over small homogeneous and large heterogeneous populations. Yashima and Sasaki in 2013 extended the MC epidemic model from over a population to a particular commute network model, for studying the epidemic spread of an infectious disease within a metropolitan area — Tokoyo train station. MC is used to simulate the spread of infectious disease by considering the commuters flow dynamics, the population sizes and other factors, the proceeding size of the epidemic and the timing of the epidemic peak. It is claimed that the MC model is able to serve as a pre-warning system forecasting the incoming spread of infection prior to its actual arrival. Narrowing from the MC model which can capture the temporal–spatial dynamics of the epidemic spread, a more specific MC model is constructed by Fitzgerald et al. [7] in 2017 for simulating queuing behaviour of an emergency department. The model incorporates queuing theory and buffer occupancy which mimic the demand and nursing resource in the emergency department respectively. It was found that adding a separate fast track helps relieving the burden on handling of patient and cutting down the overall median wait times during an emergency virus outbreak and the operation hours are at peak. Mielczarek and Zabawa [8] adopted a similar MC model to investigate how erratic the population is, hence the changes in the number of infested patients affect the fluctuations in emergency medical services, assuming there are epidemiological changes such as call-for-services, urgent admission to hospital and ICU usages. Based on some empirical data obtained from EMS centre at Lower Silesia Region in Poland, the EMS events and changes in demographic information are simulated as random variables. Due to the randomness of the changes (in population sizes as people migrate out, and infested cases increase) in both demand and supply of an EMS, the less-structured model cannot be easily examined by deterministic analytic means. However, MC model allows decision makers to predict by studying the probabilities of possible outcomes on how the changes impact the effectiveness of the Polish EMS system. There are similar works which tap on the stochastic nature of MC model for finding the most effective escape route during emergency evacuation [9] and modelling emergency responses [10].

Motivation and contributions

Overall, the above-mentioned related works have several features in common: their studies are centred on using a probabilistic approach to model complex real-life phenomena, where a deterministic model may fall short in precisely finding the right parameters to cater for every detail. The MC model is parsimonious that means the model can achieve a satisfactory level of explanation or insights by requiring as few predictor variables as possible. The model which uses minimum predictor variables and offers good explanation is selected by some goodness of fit as BIC model selection criterion. The input or predictor variables are often dynamic in nature whose values change over some spatial–temporal distribution. Finally, the situation in question, which is simulated by MC model, is not only complex but a prior in nature. Just like the new COVID-19 pandemic, nobody can tell when or whether it will end in the near future, as it depends on too many dynamic variables. While the challenges of establishing an effective MC model is acknowledged for modelling a completely new epidemic behaviour, the model reported in [10] inspires us to design the MC model by decomposing it into several sub-problems. Therefore, we proposed a new MC model called composite MC or CMC in short which accepts predictor variables from multi-prong data sources that have either correlations or some kind of dependencies from one another. This helps provide the decision maker a fuller view hence richer information that contribute to an optimal decision (so called the best out of the worst decision). The challenge here is to ensure that the input variables though they may come from random distribution, their underlying inference patterns must contribute to the final outcome in question. The methodology is focused on finding the best possible algorithmic approaches for generating stochastic future outcomes, when a novel epidemic just emerged with very limited available information on hand.

In computer science, this represents a typical problem of machine learning over incomplete/limited data in early epidemic. Considering multi-prong data sources widen the spectrum of possibly related input data, thereby enhancing the performance of Monte Carlo simulation. However, naive MC by default does not have any function to decide on the importance of input variables. It is known that what matters for the underlying inference engine of MC is the historical data distribution which tells none or little information about the input variables prior to the running of MC simulation. To this end, we propose a pre-processor, in the form of optimized neural network namely BFGS-Polynomial Neural Network is used at the front of MC simulator. BFGS-PNN serves as both a filter for selecting important variables and forecaster which generates future time-series as parts of the input variables to the MC model. Traditionally all the input variables to the MC are distributions that are drawn from the past data which are usually random, uniform or some customized distribution of sophisticated shape. In our proposed model here, a hybrid input that is composed of both deterministic type and non-deterministic type of variables. Deterministic variables come from the forecasted time-series which are the outputs of the BFGS-PNN. Non-deterministic variables are the usual random samples that are drawn from data distributions. In the case of COVID-19, the future forecasts of the time-series are the predictions of the number of confirmed infection cases and the number of cured cases. Observing from the historical records, nevertheless, these two variables display very erratic trends, one of them contains extreme outliers. They are difficult to be closely modelled by any probability density function; intuitively imposing any standard data distribution will not be helpful to delivering accurate outcomes from the MC model. Therefore in our proposal, it is needed to use a polynomial style of self-evolving neural network that was found to be one of the most suited machine learning algorithm in our prior work [11], to render a most likely future curve that is unique to that particular data variable.

The composition of the multiple data sources is of those relevant to the development (rise-and-decline) of the COVID-19 epidemic. Specifically, a case of how much daily monetary budget that is required to struggle against the infection spread is to be modelled by MC. The data sources of these factors are publicly available from the Chinese government websites. More details follow in Section 2 below. The rationale behind using a composite model is that what appears to be an important figure, e.g. the number of suspected cases are directly and indirectly related to a number of sub-problems which of each carries different levels of uncertainty: how a person gets infested, depends on 1) the intensity of travel (within a community, suburb, inter-city, or oversea) 2) preventive measures 3) trail tracking of the suspected and quarantining them 4) medical resources (isolation beds) available, and 5) eventual cured or dead. Some of these data sources are in opposing influences to one another. For example, the tracking and quarantine measures gets tighten up, the number of infested drops, and vice-versa. In theory, more relevant data are available, the better the performance and more accurate outcomes of the MC can provide. MC plays an important role here as the simulation is founded on probabilistic basis, the situation and its factors are nothing but uncertainty. Given the available data is scare as the epidemic is new, any deterministic model is prone to high-errors under such high uncertainty about the future.

The contribution of this work has been twofold. Firstly, a composite MC model, called CMCM is proposed which advocates using non-deterministic data distributions along with future predictions from a deterministic model. The deterministic model in use should be one that is selected from a collection of machine learning models that is capable to minimize the prediction error with its model parameters appropriately optimized. The advantage of using both fits into the MC model is the flexibility that embraces some input variables which are solely comprised of historical data, e.g. trends of people infested. And those that underlying elements which contribute to the high uncertainty, e.g. the chances of people gather, are best represented in probabilistic distribution as non-deterministic variables to the MC model. By this approach, a better-quality MC model can be established, the outcomes from the MC model become more trustworthy. Secondly, the sensitivity chart obtained from the MC simulation is used as corrective feedback to rules that are generated from a fuzzy rule induction (FRI) model. It is known that FRI outputs decision rules with probabilities/certainty for each individual rule. A rule consists of a series of testing nodes without any priority weights. By referencing the feedbacks from the sensitivity chart, decision makers can relate the priority of the variables which are the tests along each sequence of decision rules. Combining the twofold advantages, even under the conditions of high uncertainty, decision makers are benefited with a better-quality MC model which embraces considerations of composite input variables, and fuzzy decision rules with tests ranked in priority. This tool offers a comprehensive decision support at its best effort under high uncertainty.

The remaining paper is structured as follow. Section 2 describes the proposed methodology called GROOMS+CMCM, followed by introduction of two key soft computing algorithms — BFGS-PNN and FRI which is adopted for forecasting some particular future trends as inputs to the MC model and generating fuzzy decision rules respectively. Section 3 presents some preliminary results from the proposed model followed by discussion. Section 4 concludes this paper.

A novel methodology

MC has been applied for estimating epidemic characteristics by researchers over the year, because the nature of epidemic itself and its influences are full of uncertainty. An application that is relatively less looked at but important is the direct cost of fighting the virus. The direct cost is critical to keep the virus at bay when it is still early before becoming a pandemic. But often it is hard to estimate during the early days because of too many unknown factors. Jiang et al. [12] have modelled the shape of a typical epidemic concluding that the curve is almost exponential; it took only less than a week from the first case growing to its peak. If appropriate and urgent preventive measure was applied early to have it stopped in time, the virus would probably not escalate into an epidemic then pandemic. Ironically, during the first few days (some call it the golden critical hours), most of the time within this critical window was spent on observation, study, even debating for funding allocation and strategies to apply. If a effective simulation tool such as the one that is proposed here, decision makers can be better informed the costs involved and the corresponding uncertainty and risks. Therefore, the methodology would have to be designed in mind that the available data is limited, each functional component of the methodology in the form of soft computing model should be made as accurate as possible. Being able to work with limited data, flexible in simulating input variables (hybrid deterministic and its counterpart), and informative outcomes coupled with fuzzy rules and risks, would be useful for experts making sound decision at the critical time. Our novel methodology is based on Group of Optimized and Multisource Selection, (GROOMS) methodology [11] which is made for choosing a machine learning method which has the highest level of accuracy. GROOMS as a standalone optimizing process is in aid of assuring the deterministic model that is to be used as input variable for the subsequent MC simulation to have the most accurate data source input. By default, MC model at its naive form accepts only input variable from a limited range of standard data distributions (Uniform, Normal, Bernoulli, Pareto, etc.); best fitting curve technique is applied should the historical data shape falls out of the common data distribution types. However, this limitation is expanded in our composite MC model, so-called CMCM in such a way that all relevant data sources are embraced, both direct and indirect types. An enhanced version of neural network is used to firstly capture the non-linearity (often with high irregularity and lack of apparent trends and seasonality) of the historical data. Out of the full spectrum of data sources, direct and indirect, the selected data sources through feature selection by correlation, that are filtered by the neural network, whose data distributions would be taken as input variables to the MC model. The combined methodology, GROOMS+CMCM is shown in Fig. 1.

Fig. 1

GROOMS+CMCM methodology.

Broyden–Fletcher–Goldfarb–Shanno optimized polynomial neural network (BFGS-PNN)

According to the methodology, a machine learning algorithm candidate called BGFS-PNN which is basically PNN as selected as the winning algorithm in [11] enhanced with further parameter optimization function. The given time-series data fluctuated more than the same that were collected earlier. As a data pre-processor, BFGS-PNN has two functions. Firstly, for the non-deterministic data, using a ClassifierBasedFilter function in a wrapper approach, salient features could be found in feature selection. The selected salient features are those very relevant to the forecast target in the MC. In this case, it is composite MC model or CMCM as the simulation engine intakes multiple data sources from types of deterministic and non-deterministic. The second function is to forecast a future time-series as a type of deterministic input variable for the CMCM. The formulation of BFGS-PNN is shown as follow. The naïve version of PNN is identical to the one reported in [11]. BFGS-PNN uses BFGS (Broyden–Fletcher–Goldfarb–Shanno) algorithm to optimize the parameters and network structure size in an iterative manner using hill-climbing technique. BFGS theory is basically for solving non-linear optimization problem iteratively by finding a stationary equilibrium through Quasi-Newton method [13] and secant method [14]. Let PNN [15] take the form of Kolmogorov–Gabor polynomial as a functional series in Eq. (1).

The polynomial is capable to take form of any function which is generalized as (). The induction process is a matter of finding all the values for the coefficient vector . As the process iterates, the variables from arrive in sequence fitting into the polynomial via regression and minimizing the error. The complexity grows incrementally by trying to add a neuron at a time while the forecasting error is being monitored [16]. When the number of neurons reach a pre-set level, the hidden layer increases. This continues until there is no further performance gain observed, the growth of the polynomial stops and taken as the final equation for the PNN. It is noted that the increment of the network growth is linear.

However, for the case of BFGS-PNN, the expansion of the polynomial is non-linear and heuristic. The optimal state is achieved by unconstrainted hill-climbing method guided by Quasi-Network and secant methods. Let an error function be e(p) where p is a vector of real numbers be a vector of network structure information and parameters, i.e. neurons and layers in an ordered set. At the start, when , is initialized by randomly chosen states. Let search direction be at iteration i where time . Let be Hessian which is a square matrix of 2nd-order partial derivatives of function e. where i is the current iteration, H improves becomes a better estimate as the process iterates. e() is the gradient of the error function which needs to be minimized at by following the quasi-Newton search pattern using a gradient function similar to Eq. (2). It seeks to search for the next state of parameters values +1 by optimizing e(+S*s ) where the scalar S must be greater than 0. Eq. (2) then would have to obey the quasi-Newton condition upon solving the approximation of the Hessian Hi, in a way of Eq. (3).

By the secant equation, we let the Hessian matrix take the form as in Eq. (4).

The updating function for Hessian matrix is defined as Eq. (5) following the secant method. The equation tries to impose the condition and symmetry such that . Let , we can obtain two sub-equations in Eq. (6) as coefficients to Eq. (5). Substituting Eq. (6) Into Eq. (5), we obtain the updating function for Hessian matrix .

By applying Sherman–Morrison formula [17] to Eq. (6), we get Eq. (8) Which is the inverse of Hessian H matrix. Expanding Eq. (8) to Eq. (9) we obtain an equation that can be calculated quickly without needing any buffer space for fast optimization which aims at minimizing e(*).

Fuzzy Rule induction

By our GROOMS+CMCM methodology, raw data from multiple sources are filtered, condensed, converted into insights of future behaviours in several forms. Traditionally in MC simulation, probability density functions as simulated outcomes are generated, so is sensitivity chart which ranks how each factor in relevance to the predicted outcome. Fuzzy Rule Induction (FRI) plays a role in the methodology by inferring a rule-based model which supplies a series of conditional tests that lead to some consequences based on the same data that were loaded into the MC engine. FRI serves the threefold purpose of easy to use, neural and scalable. Firstly, the decision rules are interpretable by users. They can complement the probability density functions which show a macro view of the situation. FRI helps give another perspective in causality assisting decision makers to investigate the logics of cause-and-effect. Furthermore, different from other decision rule models, FRI allows some fuzzy relaxation in bracketing the upper and lower bounds, thereby a decision can be made based on the min–max values pertaining to each conditional test (attribute in the data). The FRI rules are formatted as branching-logic, which is also known as predicate-logic that preserves the crudest form of knowledge representation. Predicate logic has an IF-TEST[Min-Max]-THEN-VERDICT basic structure and a propensity of how often it exists in the dataset. The number of different groups of FRI rules depend on how many different labels in the predicted class. The second advantage is that the FRI rules are objective and free from human bias as they are derived homogenously from the data. Therefore, they are suitable ingredient for scientifically devising policy and strategy for epidemic control. Thirdly FRI rules can scale up or down not only in quantity, but also in cardinality. A rule can consist of as many tests as the attributes of the data are available. In other words, as a composite MC system, new source of data could be chipped in as per when it becomes necessary or recently available; the attributes of the new data can add on to the conditional tests of the FRI rules.

One drawback about FRI is the lack of indicators for each specific conditional test (or attribute). By the current formulation of FRI, the likelihood of the occurrence of the rule is assigned to the rule as a whole. Little is known about how each conditional test on the attribute is relatively contributing to the outcome that is specified in the rule. In the light of this shortcoming, our proposed methodology suggests that the scores from the sensitivity charts with respect to the relations between the attributes and the outcome, could be used at the rule by simple majority voting.

Rules are generated as a by-product of classification in data mining. The process is through fuzzification of the data ranges and the confidence factors in their effects in classification are taken as indicators. Let a rule be a series of components constraining the attributes (with outcome = Y) in the classification model building, so that they can remain valid even though the values are fuzzified. So a rule can be expressed in the predicate rule format such that each , where Y is an membership whose labels are mapped to the class labels, a membership takes a range of positive where = (-, ]. If a rule embraces a conditional test where ( ), = [, ) if it has a test (ai ). Let = [, ] when it is comprised of two tests are combined. Replacing the absolute memberships by fuzzy memberships which are defined as fuzzy sets governed by trapezoidal membership function [18]. A fuzzy trapezoidal membership normally would have to be structured by four parameters and expressed as = (). The functional logic is shown in Eq. (10): where and are the lower and upper bounds of the argument which will map to a fuzzy membership of value 1. Similarly, the supports of the lower and upper bounds are denoted by and .

The fuzzy rules are built on the decision rules which are generated from standard decision tree algorithm, such as direct rule-based classifier equipped with incremental reduced error pruning via greedy search [19]. Given the rule sets generated, the task here is to find the most suitable fuzzy extension for each rule. The task can be seen as replacing the current memberships of the rules by their corresponding fuzzy memberships, which is not too computationally difficult as long as the rule structures and the elements the same. In order to fuzzify a membership, the following formula is applied over the antecedent ( ) of rule set while considering the relevant data ; at the same time the instances from the other antecedent ( ), where j I are to be nullified. By this approach, the instances are divided into two subsets: one subset contains all the positive instances and the other subset contains all the negative instances. After that, a purity measure is used to further separate the two groups into two extremes, positive and negative subgroups, by Eq. (12). where positive instances are denoted by , and the negative instances are by , When it comes to actual operation, a certainty factor which serves as an indicator about how much the new data instances indeed belong to a subgroup, is needed to quantify the division. Followed by segregating the data into fuzzy rules by machine learning the relations from their attributes and instance values to some label class, further indicator is needed to quantify the strength of each rule. Assume we have a new test instance x, Eq. (14) computes the support of the rules of x as follow: where () is the certainty factor pertaining to the rule. The certainty factor is expressed in Eq. (15): where denotes the subsets of training instances that are labelled as. The result that is predicted by the default classifier to be one of the class labels, is the one that has the greatest value computed from the support function Eq. (14). At times, some instances x could not be classified into rule or subgroup, that happens when (x)=0 for all classes , x could be randomly assigned or temporarily placed into a special group. Otherwise, the fuzzy rules are formed, certainty and support indicators are assigned to each one of them. The indicators mean how strong the rules are with respect to the predictive power to the class label possessed by the rules.

Experiment and results

For the purpose of validating the proposed GROOMS+CMCM methodology, empirical data proceeding from The Chinese Center for Disease Control and Prevention1 (CDCP), an official Chinese government agency in Beijing, China. Since the beginning of the COVID-19 outbreak, CDCP has been releasing the data to the public and updating them daily via mainstream media Tencent and its subsidiary2 . The data come from mainly two sources: one source is known to be deterministic in nature which is harvested from CDCP in the form of time-series starting from 25 Jan 2020 to 25 Feb 2020. A snapshot of the published data is shown in Fig. 2 which are deterministic in nature as historical facts.

Fig. 2

COVID-19 data disseminated by CDCP with daily updates.

Deterministic input variables to CMC

The data collected for this experiment are only parts of the total statistics available on the website. The data required for this experiment are the numbers of people in China who have been contracted with the COVID-19 disease in the following ways: suspected of infection by displaying some symptoms, confirmed of infection by medical tests, cumulative confirmed cases, current number of confirmed cases, current number of suspected cases, current number of critically ill, cumulative number of cured cases, cumulative number of deceased cases, recovery (cured) rate % and fatality rate %. This group of time-series are subject to GROOMS for finding the most accurate machine learning technique for obtaining the forecasts as future trends under development. In this case, BFGS-PNN was found to be the winning candidate model, hence applied here for generating future trends for each of the above-mentioned records. The forecasts based on these selected data by BFGS-PNN are shown in Fig. 3. The forecasts are in turn used as deterministic input variables to the CMC model. They have relatively lowest errors in RMSE in comparison to other time-series forecasting algorithms as tested in [11]. The rationale is to use the most accurate possible forecasted inputs for achieving the most reliable simulated outcomes from MC simulation at the best effort.

Fig. 3

Forecasts of the main input variables to Monte Carlo simulation using BFGS-PNN: (a) new_daily_increase_confirmed, (b) new_daily_increase_suspected, (c) cured_rate, and (d) death_rate.

The goal of this Monte Carlo simulation experiment, which is a part of GROOMS+CMCM is to hypothetically estimate the direct cost that is needed as an urgent part of national budget planning to control the COVID-19 epidemic. Direct cost means the cost of medical resources, that includes but not limit to medicine, personnel, facilities, and other medical supplies, directly involved in providing treatments to those patients due to the COVID-19 outbreak. Of course, the grand total cost is much wider and greater than the samples experimented here. The experiment however aims at demonstrating the possibilities and flexibility of embracing both types of deterministic and non-deterministic data inputs by the composite MC methodology.

Non-deterministic input variables to CMC

The other group of data that would feed into the CMC are non-deterministic or probabilistic because they would have to bear a high level of uncertainty. They are subject to situation that changes dynamically and there is little control over the outcome. In the case of COVID-19 epidemic control, finding the cure to the virus is a good example. There is best effort put into treatment, but no certainty at all about a cure, let alone knowing when exactly a cure could be developed, tested and proven to be effective for use against the novel virus. There are other probabilistic factors which are used in this experiment as well. Selected main attributes are tabulated in Table 1. We assume a simple equation for estimating the direct cost in fighting COVID-19 using only data of quarantining and isolated medical treatment as follow. Note that the variables shown are abbreviated from the term names. E.g. = days_till_recovery. The assumptions and hypothesis are derived from past experiences about direct costs involved in quarantine and isolation during the epidemic of SARS in 2003 as published in [20], with reasonable adjustment. The reasonable adjustment is the consideration of inflation since 2003, the capacity of medical infrastructure of today’s compared to that of 2003. Another factor is the fact that during COVID-19, the Chinese government has put in unprecedent efforts in building temporary hospitals and pouring in tens of thousands of medical doctors and nurses in an attempt of curing the infested patients. This overhead cost came directly from the national funding, that should be taken account as a part of the total_daily_cost which embraces direct and indirect material costs, additional resources in the 5G communication infrastructure and a large number of hospital beds that were exclusively built for containing the COVID-19 outbreak. This effort was not seen in the same intensity and magnitude for the SARS containment. The equation for the total_daily_cost is listed in Eq. (16).

Table 1

Selected input variables used in the CMC model.

The variable total_daily_cost is a direct variable contributing to the final medical cost in this simulation scenario. However, this variable is influenced by a combination of deterministic and non-deterministic variables under high uncertainty. The values for the deterministic variables such as the cured_rate and death_rate come from the forecast of the best selected prediction algorithm. For the non-deterministic variables ppi/day, and , the following assumptions are derived from [21]. These variables are probabilistic in nature as shown in Table 1. E.g. nobody can actually tell how long an infested patient could be recovered and go home, nor how long the isolation needs to be when the patient is in critical condition. All these are bounded by some probabilities that can be expressed in statistical properties, such as min–max, mean, standard deviation and so on. So some probabilism functions are needed to describe them, and random samples from these probabilities distributions are drawn to run the simulation. The variable Days_till_death ( ) however is derived from SARS case which are not exactly the same, but the most similar in terms of uncertainty. When the situation becomes highly unpredictable and the disease is unknown in the early stage, an educated guess is made based on whatever that is found similar in the past. So the days assumed for Days_till_death is taken from the SARS case from reference [21].

It is assumed that the growth of daily medical cost ppi/day follows a normal distribution with daily increase rate. The daily increase rate is estimated from [20] to be rising as days go by, because the Chinese government has been increasingly putting in resources to stop the epidemic using national efforts. The increase is due to the daily increase number of medical staff who are flown to Wuhan from other cities, and the rise of the volume of consumable medical items as well as their inflating costs. The daily cost is anticipated to become increasingly higher as long as the battle against COVID-19 continues at full force. There are other supporting material costs and infrastructure relocation costs such as imposing curfews and economy damages. However, these other costs are not considered for the purpose of demonstrating with a simple and easy-to-understand CMC model.

Normal distribution and uniform distribution are assumed accountable for the increase and probability distributions that describe the lengths of hospital stay. When more information become available, one can consider refining them to Weibull distribution and Rayleigh distribution which can better describe the progress of the epidemic in terms of severity and dual statistical degree of freedom.

This is a very simplified approach in guessing the daily cost of the so-called medical expenses, based on only two interventions — quarantine and isolation. Nevertheless, this CMC model though simplified, is serving as an example of how Monte Carlo style of modelling can help generate probabilistic outcomes, and to demonstrate the flexibility and scalability of the modelling system. Theoretically, the CMC system can expand from considering two direct inputs (quarantining and isolation) to 20, or even 200 other direct and indirect inputs to estimate the future behaviour of the epidemic. In practical application, data that are to be considered shall be widely collected, pre-processed, examined for quality check and relevance check (via GROOMS), and then carefully loading into the CMC system for getting the outcomes.

Decision making under uncertainty

Stochastic simulation is well-suited for studying the risk factors of infectious disease outbreaks which always changes in their figures across time and geographical dispersion, thereby posing high level of uncertainty in decision making. Each model forecast by MC simulation is an abstraction of a situation under observation — in our experiment, it is the impact of the dynamics of epidemic development on the direct medical costs against COVID-19. The model forecast depicts the future tendencies in real-life situation rather than statements of future occurrence. The output of MC simulation sheds light in understanding the possibilities of outcomes anticipated.

Being a composite MC model, the ultimate performance of the simulated outcomes would be sensitive to the choice of the machine learning technique that generated the deterministic forecast as input variable to the CMC model. In light of this, a part of our experiment besides showcasing the MC outcomes by the best available technique, is to compare the levels of accuracy (or error) resulted from the wining candidate of GROOMS and a standard (default) approach. The forecasting algorithms in comparison are BFGS+PNN and Linear Regression respectively. The performance criterion is RMSE, which is consistent and unitless and RMSE=, as defined in [22]. At the same time, the total costs that are produced manually by explicitly use of spreadsheet using human knowledge are compared vis-à-vis with those of the forecast models by CMCM. The comparative performances are tabulated in Table 2. The forecasting period is 14 days. The CMCM model is implemented on Oracle Crystal Ball Release 11.1.2.4.850 (64-bit), running on a i7 CPU @ 2 GHz, 16Gb RAM and MS Windows 10 platform. 10,000 trials were set to run the simulation for each model. The reason for choosing 10,000 trials is to try obtaining a stable result. The 10k number is an educated guess. The typical number of trials as suggested by the maker of the software is 5,000. For a good safety margin, the recommended number for a result is doubled considering the comprise of the length of simulation time and the incremental differences of the results when scaling up the number of trials. It was found that by scaling over 10,000 trials and more, the results differ in less than 0.01%. The parameters that are used to run the MCMC model are shown in Appendix.

Table 2

Comparative performance of various forecast model approaches.

It is apparent that as seen that from Table 2, the RMSE of the Monte Carlo forecasting method using linear regression is more than double that of the method using BFGS+PNN (approximately 128K vs 62K). That is mainly due to the overly estimate of all the deterministic input variables by linear regression. Referring to the first diagram in Fig. 3, the variable called new_daily_increase_confirmed is non-stationary and it contains an outlier which rose up unusually high near the end. Furthermore, the other correlated variable called new_daily_increase_suspected, which is a precedent of the trends of the confirmed cases, also is non-stationary and having an upward trend near the end, though it dips eventually. By using linear regression, the outlier and the upward trends would encourage the predicted curve to continue with the upward trend, linearly, with perhaps a sharp gradient. Consequently, most of the forecast outcomes in the system have been over-forecasted. As such, using linear regression causes unstable stochastic simulation, leading to the more extreme final results, compared to the other methods. This is evident that the total_daily_cost has been largely over-forecasted and under-forecasted by manual and MC approaches, in Table 2. Nevertheless this observation was accurate only up to the time the article was written. New developments of pandemic hence the data evolve every day, the results and conclusions are anticipated to be different in future days to come.

On the other hand, when the BFGS+PNN is used, which is able to better recognize non-linear mapping between attributes and prediction class in forecasting, offers more realistic trends which in turn are loaded into the CMCM. As a result, the range between the final total_daily_cost results are narrower compared it to its peer linear regression ([LR: 12mil – 83mil] vs [BFGS+PNN: 54mil – 74mil]). The estimated direct medical cost for fighting COVID-19 for a fortnight is estimated to be about 73.6 million USD given the available data using GROOMS+CMCM.

According to the results in the form of probability distributions in Fig. 4, different options are available for the user to choose from when it comes to estimating the fortnightly budget in fighting this COVID-19 given the existing situation. Each option comes with different price tag, and at different levels of risks. In general, the higher risk that the user thinks it can be tolerated, the lower the budget it will be, and vice-versa. From the simulated possible outcomes in Fig. 4, if budget is of constrained, the user can consider bearing the risk (uncertainty of 50%) that a mean of $74mil with [min:$69mil, max:$78mil] is forecasted to be sufficient to fulfil the direct medical cost need. Likewise, if a high certainty is to be assured, for example at 80% the chance that the required budget would be met, it needs about a mean of $79mil with [min:$66mil, max:$82mil]. For a high certainty of 98%, it is forecasted that the budget will fall within a range of [min:$60mil, max:$89mil] at the mean of $80mil. As a de-facto practice, some users will take 80% certainty as Pareto Principle (80–20) decision [23] and accept the mean budget at $79mil. $79mil should be a realistic and compromising figure when compared to manual forecast without stochastic simulation, where $54mil and $84mil budgets would have been forecasted by manual approach by linear regression and neural network respectively.

Fig. 4

Probability Distributions as outcomes of CMCM at certainly levels of (a) 50%, (b) 80% and (C) 98%.

Sensitivity chart and Fuzzy rules

Sensitivity chart, by its name, displays the extents of how the output of a simulated MC model is affected by changes in some of the input variables. It is useful in risk analysis of a so-called black box model such as the CMCM used in this experiment by opening up the information about how sensitive the simulated outcome is to variations in the input variables. Since the MC output is an opaque function of multiple inputs of composite variables that were blended and simulated in a random fashion over many times, the exact relationship between the input variables and the simulated outcome will not be known except through sensitivity chart. An output of sensitivity chart from our experiment is generated and shown in Fig. 5.

Fig. 5

Sensitivity chart showing the measured sensitivity between each variable and the output.

As it can be observed from Fig. 5, the top three input variables which are most influential to the predicted output, which is the total medical cost in our experiment are: the average number of days before recovery in Day 10 and Day 12, and average cost per day for isolating a patient in day 2. The first two key variables are about how soon a patient can recover from COVID-19, near the final days. And the third most important variable is the average daily cost for isolating patients at the beginning of the forecasting period. This insight can be interpreted that an early recovery near the final days and reasonably lower medical cost at the beginning would impact the final budget to the greatest extend. Consequently, decision makers could consider that based on the results from the sensitivity analysis, putting in large or maximum efforts in treating isolating patients at the beginning, observe for a period of 10 days or so; if the medical efforts that were invested in the early days take effect, the last several days of recovery will become promising, hence leading to perhaps saving substantially a large portion of medical bill.

Sensitivity chart can be extended to what-if scenario analysis for epidemic key variable selection and modelling [24]. For example, one can modify the quantity of each of the variables, and the effects on the output will be updated instantly. However, this is beyond the scope of this paper though it is worth exploring for it is helpful to fine-tune how the variables should be managed in the effort of maximizing or minimizing the budget and impacts.

Since the effect of a group of independent variables on the predicted output is known and ranked from the chart, it could be used as an alternative to feature ranking or feature engineering in data mining. The sensitivity chart is a by-product generated by the MC after a long trial of repeating runs using different random samples from the input distributions. Fig. 6 depicts how the sensitivity chart is related to the processes in the proposed methodology. Effectively the top ranked variables could be used to infer the most influential or relevant attributes from the dataset that loads into an FRI model (described in Section 2.2) for supervised learning. One suggested approach which is fast and easy is to create a correlogram, from there one can do pairwise matching between the most sensitive variables from the non-deterministic data sources to the corresponding attributes from the deterministic dataset. Ranking scores could be mapped over too, by Boyer–Moore majority vote algorithm [25].

Fig. 6

Sensitivity chart and fuzzy rule induction as outputs from the GROOMS-CMCM methodology (*accuracy by GROOMS).

Some selected fuzzy rules generated by the methodology and filtered by the sensitivity chart correlation mapping [26] are show below. The display threshold is 0.82 which is arbitrarily chosen to display only the top six rules where half of them predict a reflection point for the struggle of controlling the COVID-19 can be attained, the other half indicate otherwise. An FRI model in a nutshell is a classification which predicts an output belonging either one of two classes. In our experiment, we set up a classification model using FRI to predict whether an inflection point of the epidemic could be reached. There is no standard definition of inflection point, though it is generally agreed that is a turning point at which the momentum of accumulation changes from one direction to another or vice-versa. That could be interpreted as either an intersection of two curves of which their trajectories begins to switch. In the context of epidemic control, an inflection point is the moment since when the rate of spreading starts to subside, thereafter the trend of the epidemic is leading to elimination or eradication.

Based on a sliding window of 3 days length, a formula for computing the inflection point based on the three main attributes of the COVID-19 data is given in Eq. (18), which is listed as follow: where , , and which can be arbitrarily set by the user. The weights reflect the importance which one considers on how the up or downward trends of confirmed cases and cured vs death rates contribute to reaching the inflection point. A dual curve chart that depicts the inflection point is shown in Fig. 7.

Fig. 7

Curve chart that approximates the search for inflection point during the epidemic.

Interestingly, when near the end of the timeline (28/1 – 20/2) that is from point 19th onwards, the two curves seem to intervene, as it has been hoping that the winning curve is rise over the losing curve. An inflection point might have reached, but the momentum of the winning is there yet. Further observation on the epidemic development is needed to confirm about the certainty of winning. Nevertheless, the top six rules that are built from classification of inflection point, and processed by feature selection via sensitivity analysis, are shown below. CF stands for confidence factor which indicates how strong the rule is. On the winning side, the rules reveal that when the variables about new confirmed cases fall below certain numbers, a win is scored, contributing towards a reflection point.

Win:

RULE 1: (yester3days-ndic = ’(-.. 3581]’) win=1 (CF = 0.96)

RULE 2: (ndic = ’(-.. 1874.4]’) win=1 (CF = 0.92)

RULE 3: (yester2days-ndic = ’(-.. 5027.5]’) win=1 (CF = 0.85)

Lose:

RULE 4: (cured_rate = ’(-.. 3.86]’) & (ndic = ’(1874.4 .. 3349.8]’) win=0 (CF = 0.8)

RULE 5: (ndic = ’(13677.6 .. )’) win=0 (CF = 0.53)

RULE 6: (yester3days-ndic = ’(10706.5 .. )’) win=0 (CF = 0.53)

The strongest rules of the two forces are Rules 1 and 4. Rule 1 shows that to win an epidemic control the down trend over consecutive three days must fall below 3581; on the other hand, the epidemic control may lead to failure, if the cured rate stays less than 3.86% and the new daily increase in confirmed cases remain high between 1874 and 3350 (round up the decimal points).

Conclusion

Since the outbreak of the novel coronavirus (COVID-19) in December 2019, the pandemic which started as a very contagious domestic epidemic gave rise to an emergency situation where critical decisions were demanded for. The virus was novel and very little information was known about the epidemic at the early stage. With incomplete information, limited data on hand, and ever changing on the epidemic development, it is extremely hard for anybody to decide using only deterministic approach which foretells precisely the future behaviour of the epidemic. The work reported in this paper is focused on solving the problem of incomplete data available for stochastic decision making. In this paper a composite Monte-Carlo model (CMCM) is proposed to be used in conjunction with GROOMS methodology [21] which finds the best performing deterministic forecasting algorithm. Coupling GROOMS+CMCM together offers the flexibility of embracing both deterministic and non-deterministic input data into the Monte Carlo simulation where random samples are drawn from the distributions of the data from the non-deterministic data sources for reliable outputs. During the early period of disease outbreaks, data are scarce and full of uncertainty. The advantage of CMC is that a range of possible outcomes are generated associated with probabilities. Subsequently sensitivity analysis, what-if analysis and other scenario planning can be done for decision support. As a part of the GROOMS+CMCM methodology, fuzzy rule induction is also proposed, which provides another dimension of insights in the form of decision rules for decision support. A case study of the recent novel coronavirus epidemic (which are also known as Wuhan coronavirus, COVID-19 or 2019-nCoV) is used as an example in demonstrating the efficacy of GROOMS+CMCM. Through the experimentation over the empirical COVID-19 data collected from the Chinese government agency, it was found that the outcomes generated from Monte Carlo simulation are superior to the traditional methods. A collection of soft computing techniques, such as BFGS+PNN, Fuzzy Rule Induction, and other supporting algorithms to GROOMS+CMCM could be able to produce qualitative results for better decision support, when used together with Monte Carlo simulation, than any of deterministic forecasters alone. Forecasting carries very high uncertainty; the forecast results are much influenced by government policy and human intervention factors. The essence is to obtain the most accurate possible results given what is available in the current situation. The main contribution of our methodology is not to pursue an utmost​ accurate forecasting model; rather, it is on about fusing deterministic and non-deterministic data series into a MC simulation for fuzzy decision making. The experiment reported in this paper represents a typical scenario of decision making for a novel disease in the early days since the outbreak. It does not mean to describe how this COVID-19 progressed, as it still evolves at the time of writing. However, it shows a possible solution in generating future outcomes from the GROOMS+CMCM when the data is limited and incomplete. Early decision making is critical in the initial stage of novel epidemic, using multi-prong data sources, from both deterministic and non-deterministic domains. GROOMS+CMCM was designed to fulfil this solution.

Radar chart of the performance comparison.

Future works

By the time this article was written, there are thousands of people struggling for lives in hospital and WHO has already declared COVID-19 as a global pandemic. Latest data are updated frequently [27]. It is anticipated that MCMC would evolve in several directions, becoming more useful, detailed and pragmatic in actual deployment. One direction is the fusion of big data analytics which is going to be a truly multi-prong, multi-dimensional, and full volume of all possibly available information to be loaded into MCMC for the most reliable simulation. Following this direction, several technical challenges are required to be met. That includes extraction of only the most relevant information from big data which are often sparse and unstructured, efficient processing cloud platform and powerful computing hardware [28], and theoretical validation of the MCMC model as well as its results [29]. There are similar prior arts [30], [31] in modelling epidemic reported in literature. An ensemble like MC model [32] is applied to modelling a group of components using different MC model, and piecing up the results together for a full view. Although their underlying techniques are different, it is academically interesting to compare and contrast the differences between them, as well as their pros and cons. Most of the past works targeted at modelling disease itself. The methodology proposed here is to enhance decision making for novel disease. The contribution in this paper is unique for its one of the few approaches using stochastic simulation using multi-prong data sources. They are important because the decision made could be directly or indirectly determining the lives of thousands of patients and the fate of the epidemic control. Other directions are pertaining to the input data which shall be refined to micro-level for accuracy. Micro-level data include those influential data that describe the epidemic events, e.g. the demography of the patients, their prognosis, and records of critical illness, if any; precise medical costs per day per hospital, and the mobility of people who travels which shall be modelled accurately, for the sake of generating highly possible outcomes.

CRediT authorship contribution statement

Simon James Fong: Conceptualization, Software. Gloria Li: Data curation. Nilanjan Dey: Investigation. Rubén González Crespo: Investigation, Validation, Writing - review & editing. Enrique Herrera-Viedma: Writing - review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Table 3

Assumed parameters values for the simulation model.

Methods	Rank	RMSE	MAD	MAPE	Theil’s U	Durbin–Watson
ARIMA(0,1,2)	3	0.377797286	0.154576397	0.075464854	0.525839954	1.992976321
ARIMA(2,1,1)	2	1.650040667	0.822540658	0.144476662	0.973748908	2.187529197
ARIMA(5,2,2)	1	539.9901609	334.6537318	0.122208376	0.707004411	2.016801928
ARIMA(9,1,4)	1	1982.877409	1181.854077	0.56389354	1.596011778	1.967498782
Damped Trend Non-Seasonal	3.25	736.3747536	365.0306825	0.240931097	0.982062458	2.033016783
Double Exponential Smoothing	3.75	746.3788246	410.31494	0.241303134	1.044145056	1.955115204
Double Moving Average	3.25	949.1562194	562.5036533	0.245757604	2.044582566	1.065237993
Single Exponential Smoothing	4.5	764.4758927	398.057743	0.258685136	1.031872233	1.977003965
Single Moving Average	4.5	828.4515473	425.6190098	0.250385153	1.429116336	1.366212845

Table 4

Performance in RMSE for each key input elements to MCMC using different classical time series forecasting methods and data mining methods.

ppi/day:	34.93
ppq/day:	25.12
min_days_for_recovery:	11
max_days_for_recovery:	26
mean_days_till_death:	35.9
std_days_till_death:	6.37
quarantine_days:	14
last_current_suspected:	7264
last_current_confirmed:	58010
last_cured_rate:	15.2
last_death_rate:	2.5
ppi_daily_increase_rate:	0.005
std_ppi/day:	9