Mean Monte Carlo Finite Difference Method for Random Sampling of a Nonlinear Epidemic System

Abstract

In this article, a numerical method integrated with statistical data simulation technique is introduced to solve a nonlinear system of ordinary differential equations with multiple random variable coefficients. The utilization of Monte Carlo simulation with central divided difference formula of finite difference (FD) method is repeated n times to simulate values of the variable coefficients as random sampling instead being limited as real values with respect to time. The mean of the n final solutions via this integrated technique, named in short as mean Monte Carlo finite difference (MMCFD) method, represents the final solution of the system. This method is proposed for the first time to calculate the numerical solution obtained for each subpopulation as a vector distribution. The numerical outputs are tabulated, graphed, and compared with previous statistical estimations for 2013, 2015, and 2030, respectively. The solutions of FD and MMCFD are found to be in good agreement with small standard deviation of the means, and small measure of difference. The new MMCFD method is useful to predict intervals of random distributions for the numerical solutions of this epidemiology model with better approximation and agreement between existing statistical estimations and FD numerical solutions.

Keywords

ordinary differential equation prediction interval finite difference method Monte Carlo method epidemiology random sampling

Differential equation arises in formulation of many real life and industrial-based dynamical relations that describe growth, transformation, and change of specific dependent and independent mechanisms that may affect the whole system of inquisition. Generally, differential equations are categorized as ordinary differential equations (ODEs) or partial differential equations (PDEs). ODEs are differential equations with only one independent variable that reflected in the denominator expression, while PDEs have two or more independent variables in the system. There are three types of solutions for differential equations, which are exact, analytical, and numerical solutions. Many authors have conducted study on differential equations using vast choices of method. Since exact solutions do not exist for many complex models, closed form and analytical solutions are preferred as second-choice alternative depending on the versatility, practicality, and mathematical severity of the problems. Mohyud-Din, Noor, and Noor (2009) reviewed the applications of the analytical methods of He’s and modified variational iterations, homotopy perturbation, parameter expansion, and exp-function to solve a wide class of nonlinear initial and boundary value problems. Further discussion on such methods is also studied in Mohyud-Din and Ali (2016), Mohyud-Din, Yildirim, and Gulkanat (2012), and Ul Hassan and Mohyud-Din (2016). On the other hand, numerical methods are favored by numerical scientists and engineers due to simpler algorithm and computer-handling operations on complex solutions, less time consumption, and when the exact or analytical solutions of certain models are not yet found.

Finite difference (FD) method is one of established numerical methods to solve differential equations approximately. This method is widely applied in practical mathematical and physical modeling (Karris 2007), numerical analysis/simulation, and discrete domain problems throughout many fields of science, engineering, social science, medicine, and economics. Witnessing continuous demand and progress in these fields, many authors employed FD in their works including on ODEs. Rao and Chakravarthy (2013) presented FD method for singularly perturbed linear second-order differential-difference equations of convection–diffusion type with a shift parameter when the shift is smaller or greater than a perturbation parameter. Moghaddam and Mostaghim (2013) used a new numerical method that depends on FD to provide feasible solution for boundary problem of delayed differential equations of fractional order. In 2014, linear second-order differential-difference equation of convection–diffusion type with singularly perturbed boundary value problem is discussed, where a tridiagonal FD method is introduced with fitting parameter from the theory of singular perturbation (Rao and Chakravarthy 2014). Recent works on the FD method can be found in Liu and Chen (2014), Su, Li, and Ding (2013), and Zhao (2013).

Monte Carlo (MC) method is an eminent statistical modeling method in applied mathematics to develop models of phenomena with unknown probabilistic property. Because of its simplicity and good statistical results, the MC method has become a widespread technique used for doubt quantity. Due to increasing demand of statistical modeling and fast development of technology, MC method is employed with multiprocessor computing system in order to simulate many independent statistical experiments at a time. This progress enables difficult problems of probabilistic nature to be solved using the MC method by direct simulation of the natural probability model as well as to solve deterministic problems that are represented as constructed probability processes. These processes are modeled, and the numerical solutions are framed in the form of statistical estimators of deterministic equations with either random coefficients or random boundary conditions or right sides random (Rubinstein and Kroese 2007). The MC method can also be applied to solve problems in finance, science, engineering, economics, and actuarial mathematics. It can be used to generate random data in behavioral studies as well as to estimate numerical quantities via repeated sampling. Moreover, mathematical algorithms can be randomized using computers to carry out statistical sampling in optimization problems via the MC method.

There are many authors who organized research in parameters estimation of ODEs with MC. Spigler and Zerbetto (2013) studied the oscillatory behavior of tension leg platform offshore structure using MC method, while De Souza, Colaço, and Leiroz (2014) provided MC solution of ODEs from a first law analysis and calculated the stochastic moments. Licea, Villafuerte, and Chen-Charpentier (2013) dealt with Riccati equation that has random coefficients and compared the solutions with MC method. Buezas, Rosales, and Sampaio (2013) solved a nonlinear ODE using finite element method with MC simulation technique. Recently, Leander, Lundh, and Jirstrand (2014) presented an article related to estimated parameters of ODEs that have experimental data for discrete time. It is demonstrated that transforming an ODE to a stochastic differential equation as an objective function eliminates problem of local minima. Rohaninejad and Zarghami (2012) predicted the behavior of impounded dams by employing MC with FD method, where the statistical test Kolmogorov–Smirnov is used to predict the dam behavior. Other applications of MC method are given in Abdulle and Blumenthal (2013), Dereich and Heidenreich (2011), and Kovtanyuk, Botkin, and Hoffmann (2012).

Obesity and overweight are excessive fats in one’s body measured using a simple formula of body mass index (BMI). These abnormal body weights have been medically proven as major risks to modern population with known linkage to many death cases around the world due to cardiovascular diseases, such as stroke and heart disease, diabetes, musculoskeletal disorders especially osteoarthritis, and also due to colon, endometrial, and breast cancers. Based on statistics from the World Health Organization, the world obesity rate has increased by 200 percent since 1980 with 42-million children under five years old are obese or overweight in 2013 alone. Mainly caused by imbalance of daily nutritional intake and lack of physical activities to burn extra calories taken, the excess weight continues to become a serious health problem, especially in the low- and middle-income countries exasperated by insufficient prevention strategies, law policies, and power enforcement by the corresponded governments. As much as the rate of excess weight is rising at worrying stage, the research and studies devoted to this issue also keep rising.

Santonja et al. (2012) studied the effects of a public campaign to reduce excess weight in a Spanish region within a community of Valencia that is caused by unhealthy lifestyles. A mathematical model is presented to track on the progression of excess weight loss influenced by prevention strategy of public health campaigns. The campaigns focus on overweight and/or obese individuals in order to reduce their weight. In Snyder (2007), the author reached to a conclusion that nutritional campaigns should be able to change feeding behavior. Christakis and Fowler (2007) investigated similar issue by assuming the excess weight (obesity and overweight) as a socially transmitted epidemic disease. Further studies on excess weight from diverse perspectives are conducted in Livingston et al. (2001); Navarro-Barrientos, Rivera, and Collins (2011); and Thomas et al. (2014).

In the present study, FD method integrated with MC technique is considered to simulate coefficients of time t (week) of a nonlinear system of equations representing three Valencia community subpopulations based on BMI (Santonja et al. 2012). The system consists of three nonlinear ODEs with multiple coefficients/parameters that are described as random variables with probability distributions. The integrated method we named as mean Monte Carlo finite difference (MMCFD) method is used to simulate the parameters of the model (Santonja et al. 2012) and secondly to calculate the prediction interval for the numerical solution of every dependent variable of the system when the real-valued coefficients are provided previously. Because of random sampling and rounding errors, the solution must be supported by prediction interval based on percentile of variables distribution. In particular, we are using the 5th and 95th percentiles in the present work. This interval of numerical results is later compared with the interval of statistical predictions (Santonja et al. 2012) derived by Latin hypercube sampling (LHS) technique. In order to obtain empirical prediction intervals, this article leads to set up interval by using the percentiles from the simulated values. This predicted interval represents empirical prediction intervals. The MMCFD solutions are also compared with the results generated by the FD method via central divided difference formula. The MMCFD method shows better accuracy and convergence to the given predicted results than the FD method. Thus, it is expected that the MMCFD method has a potential in future estimations with probabilistic model.

Mathematical Model

Consider the epidemiological model of a group of people in the Valencia community of Spain (Santonja et al. 2012) on the effects of public health campaigns toward the community weight loss. The selected group of adults at the age of 23 years old is divided into three subpopulations according to their BMI, which is calculated as BMI = weight/height²:

N—subpopulation of individuals with normal weight (BMI < 25),

S—subpopulation of individuals with overweight (25 ≤ BMI < 30), and

O—subpopulation of individuals with obesity (BMI ≥ 30).

The transition of individual weight is observed at the age of 24 years old in all subpopulations of N, S, and O, which is a year prior to the public health campaign organized in 2000. The initial subpopulations of individuals at the age of 23 years old with normal weight, overweight, and obesity are denoted by $N (t = 0) = N_{0} = 0.704$ , $S (t = 0) = S_{0} = 0.250$ , and $O (t = 0) = O_{0} = 0.046$ (Santonja et al. 2012), respectively, where t is a time unit in week with t = 0 is the last week of the year 2000. The average subpopulation proportions are represented by the terms µN, µS, and µO, respectively. It is further assumed that when the individuals reach 65 years old, they will leave the system in all subpopulations. Due to unhealthy social lifestyle, the rate of transformation from normal weight subpopulations to overweight and obese subpopulations is nonlinearly modeled by βN (S + O), while the rate of transformation of individuals from overweight to obese subpopulation is a linear term modeled by γS. Moreover, individuals in S transit to N at the rate of $ρ = ρ_{1} \times ρ_{2} \times p$ , while the individuals in O transit to S at the rate of $∊ = ∊_{1} \times ∊_{2} \times p$ . The size of the S and O subpopulations are modeled by ρS and O, respectively. Hence, the transition model is presented as a nonlinear system of ODEs with respect to independent time variable, t (week), as follows (Santonja et al. 2012):

N^{'} (t) = μ N_{0} - μ N (t) - β N (t) [S (t) + O (t)] + [p_{1} \times ρ_{2} \times ρ] S (t),

S^{'} (t) = μ S_{0} + β N (t) [S (t) + O (t)] - [μ + γ + ρ_{1} \times ρ_{2} \times ρ] S (t) + [∊_{1} \times ∊_{2} \times ρ] O (t),

O^{'} (t) = μ O_{0} + γ S (t) - [μ + ∊_{1} \times ∊_{2} \times ρ] O (t) .

Other parameters of the model (Santonja et al. 2012) are listed in Table 1 with further estimation details. This model is a system of ODEs with variable solutions depending on time t. The model parameters are random variables with estimated values. This is where the MMCFD method can be very useful to solve the model numerically.

Table 1.

Model Parameters Provided by Santonja et al. (2012).

Parameters	Remarks	Parameters	Remarks
N ₀ = 0.704	Normal weight proportion of 23-year old in 2000	p = 1/7 (per week)	Average overweight/obese individual needs to reduce a mean of 7 kg weight to transit to normal weight/overweight
S ₀ = 0.25	Overweight proportion of 23-year old in 2000	ρ	Rate that overweight individuals turn into normal weight
O ₀ = 0.046	Obese proportion of 23-year old in 2000	ρ₁ = 0.020	Two-percent proportion of overweight population who do physical activity to reduce weight
μ = 0.0004578	Average stay time in system of 24 to 65 years old adults; μ = 1/2,184 week⁻ ¹ spent time of individuals in system (42 years = 2,184 weeks)	ρ₂ = 0.042	4.2 Percentage proportion of overweight population who improve diet to reduce weight
β = 0.001121	Transmission rate that depends on unhealthy lifestyles	∊	Rate at which obese individual become overweight.
γ = 0.0003226	Rate at which 24 to 65 years old overweight adults become obese due to unhealthy lifestyles	∊₁ = 0.004	0.4 Percentage proportion of obese population who do physical activity in order to reduce weight
1cp	Average time an individual needs to return to S(t) from O(t) or N(t) from S(t)	∊₂ = 0.004	2.4 Percentage proportion of obese population who improve diet to reduce weight

Solution Approach

The general form of a nonlinear initial value problem is written as:

y^{'} = f (t, y), a \leq t \leq b, y (a) = a_{0} .

We divide [a, b] into T subintervals such as t_i = a + ih for i = 0, 1, 2,…, T, where T is the maximum number of iterations and h is a step size. The step size h = 1 (one week) is chosen here since we are analyzing the actual social epidemic model of body weight change which is more suitable to be observed by weekly at the least instead of by daily. Therefore, the equation (4) becomes:

y^{'} (t_{i}) = f (t_{i}, y (t_{i})) .

By using central difference formula of FD method, the equation (5) becomes:

\frac{d y}{d t} \approx \frac{y_{i + 1} - y_{i - 1}}{2 h}, for each i = 0, ..., T .

Upon substitution of the equation (6) into the nonlinear system of ODEs (1) to (3) gives:

\frac{N_{i + 1} - N_{i - 1}}{2 h} = μ N_{0} - μ N_{i} - β N_{i} (S_{i} + O_{i}) + (ρ_{1} \times ρ_{2} \times ρ) S_{i},

\frac{S_{i + 1} - S_{i - 1}}{2 h} = μ S_{0} + β N_{i} (S_{i} + O_{i}) - (μ + γ + ρ_{1} \times ρ_{2} \times ρ) S_{i},

\frac{O_{i + 1} - O_{i - 1}}{2 h} = μ O_{0} + γ S_{i} - (μ + ∊_{1} \times ∊_{2} \times ρ) O_{i} .

All model parameters (Santonja et al. 2012) are simulated by MC method prior to numerical solutions of the system (7) to (9) by using FD method in MMCFD approach corresponding to specific time interval (Table 2). This process is repeated for 100; 500; 1,000; 3,000; and 5,000 times. As a result, the system has 100 to 5,000 final solutions that are used to calculate their mean. The MMCFD solutions represent the system solution. These mean values correspond to the mean of the last iteration of MMCFD solutions in each repeated simulation cum the mean of the last component of final solutions vector (the end values in the ordered random sample for N(t), S(t), and O(t)) are listed in Table 2. The flow chart of MMCFD procedure is further presented in Figure 1.

Table 2.

The Mean Monte Carlo Finite Difference (MMCFD) Final Solutions and Measure of Difference for All Subpopulations.

Subpopulation	Method	2013		2015		2030
Subpopulation	Method	Solutions	\|E\|	Solutions	\|E\|	Solutions	\|E\|
N(t)	From Santonja et al. (2012)	.4350	—	.4244	—	.3658	—
	FD	.573849	.138849	.555952	.131552	.445840	.080040
	MMCFD
	100 Simulation	.564445	.12944 5	.543390	.118990	.402395	.036595
	500 Simulation	.562819	.127819	.541705	.117305	.401929	.036129
	1,000 Simulation	.563993	.128993	.543033	.118633	.403929	.038129
	3,000 Simulation	.564385	.129385	.543536	.119136	.405434	.039634
	5,000 Simulation	.564039	.129039	.543108	.118708	.404293	.038493
S(t)	From Santonja et al. (2012)	.3970	—	.4004	—	.4132	—
	FD	.325999	.071001	.335716	.064684	.388068	.025132
	MMCFD
	100 Simulation	.329961	.067039	.341145	.059255	.404574	.008626
	500 Simulation	.331786	.065214	.343118	.057282	.406629	.006571
	1,000 Simulation	.330428	.066572	.341593	.058807	.404404	.008796
	3,000 Simulation	.329959	.067041	.341030	.059370	.403192	.010008
	5,000 Simulation	.330336	.066664	.341476	.058924	.404109	.009091
O(t)	From Santonja et al. (2012)	.1680	—	.1752	—	.2210	—
	FD	.100152	.067848	.108332	.066868	.166092	.054908
	MMCFD
	100 Simulation	.104405	.063595	.114111	.061089	.190559	.030441
	500 Simulation	.104078	.063922	.113681	.061519	.188752	.032248
	1,000 Simulation	.104245	.063755	.113858	.061342	.188956	.032044
	3,000 Simulation	.104368	.063632	.113973	.061227	.188750	.032250
	5,000 Simulation	.104218	.063782	.113818	.061382	.188732	.032268

Figure 1.

Flow chart of mean Monte Carlo finite difference method.

We compute the time t as a unit in week by using Wolfram Alpha program, where the number of weeks from the beginning of 2001 to the end of 2013 is calculated as 678 weeks, from the beginning of 2001 to the end of 2015 is 782 weeks, and from the beginning of 2001 to the end of 2030 is 1,565 weeks. The number of weeks in each interval is the same number of iterations. Random predictions in 2013 and over the next few years in 2015 and 2030 can be performed by assuming all model parameters to have uniform probability distributions on given intervals (Santonja et al. 2012) since these parameters are random variables. In this case, the outcome results are called as “numerical-simulated solutions” by means of MMCFD method. Since we have previous results for this model that are estimated by a statistical method with estimated parameters (Santonja et al. 2012), comparison between the present MMCFD final solutions with the predicted results (Santonja et al. 2012) is possible. The numerical-simulated solutions of MMCFD are appeared as follows:

s o l_{n, j} (t) = {(\begin{matrix} s o l_{11} & \dots & s o l_{1 n} \\ ⋮ & ⋱ & ⋮ \\ s o l_{T 1} & \dots & s o l_{T n} \end{matrix})}_{T \times n}, where j = N, S, O,

such that each variable N(t), S(t), and O(t) is presented in the above matrix. Then, the solution of the system is:

s o l_{n} (t) = (s o l_{n, N} (t), s o l_{n, S} (t), s o l_{n, O} (t)),

when T is the maximum number of iterations as well as the number of weeks in each considered intervals (678 weeks in 2013, 782 weeks in 2015, and 1,565 weeks in 2030), and n is the number of random sample space as well as the number of simulation.

The solutions for each variable N(t), S(t), and O(t) in the model using n simulations (taking from the last iteration of equation 11) are represented as the following vector:

s o l_{n, j} (T) = {[s o l_{n} (T, i)]}_{T \times n}, where i = 1, 2, \dots, n and j = N, S, O .

Therefore, the solution of the model is:

s o l_{s y s t e m} (T) = (s o l_{n, N} (T), s o l_{n, S} (T), s o l_{n, O} (T)) .

Hence, the mean of the solutions for each system variable μ_sol,j is displayed as:

μ_{s o l, j} = \sum_{i = 1}^{n} \frac{s o l_{n} (T, i)}{n}, such that j = N, S, O .

where the MMCFD solutions of the system, μ_sol, is:

μ_{s o l} = (μ_{s o l, N}, μ_{s o l, S}, μ_{s o l, O}) .

All MMCFD computations are done in MATLAB 13 environment. Some graphical results (for 5,000 simulations) are sketched using MagicPlot (Version 2.7.2) software, while the box plots are drawn by using S-plus 8.1 software.

Statistical Indicators

This study considers certain statistical indicators such as measure of difference |E|, standard deviation of MC (σ_n), and percentile of distribution in order to compare between MMCFD results and the predicted results (Santonja et al. 2012). For this purpose, we define a measure of difference |E|:

| E | = | x_{p} - x_{s o l} |,

where x_sol is a numerical value (the last iteration value of FD or the MMCFD solutions, μ_sol,j, as in equation 14) and x_p is a predicted value from Santonja et al. (2012). Standard deviation σ is calculated as a special measuring of error for the MMCFD method. The expression $σ / \sqrt{n}$ reduces the error with the square root of the sample size n, approximating to the standard deviation of the means, σ_n, which has the following formula:

σ_{n} \approx \frac{σ}{\sqrt{n}},

where n represents the number of simulations in the present computation as well as the sample size of the model random variables. Since value of the standard deviation σ of the population is not given, the variance of σ_n defined as:

σ_{n}^{2} = 〈 I^{2} 〉 - {〈 I 〉}^{2}

is used in the current article with

〈 I 〉 = \frac{1}{n} \sum_{i = 1}^{n} {(s o l_{n, j} (T))}_{i}

and

〈 I^{2} 〉 = \frac{1}{n} \sum_{i = 1}^{n} {(s o l_{n, j} (T))}_{i}^{2},

such that j = N, S, O, where sol_n,j(T) is mentioned in equation (13). In general, the best result is obtained when the standard deviation of the means σ_n is the smallest (Jacoboni and Lugli 1989).

Another statistical indicator to describe data (solution) is prediction interval. Prediction interval is an interval contains upper and lower bounds for predicted values of the distribution for each subpopulation of the model. It can be obtained by using the pth percentile to give upper and lower limits. A percentile is a value inside a distribution that can divide ordered predicted values into two or more parts by drawing lines between these values. These predicted values belong to the vector distribution of random variable N(t), S(t), and O(t). As a consequence, the pth percentile of the predicted values is inside a population. This percentile value is equal to or less than the number one required to calculate it, when 0 < p < 100. Suppose that n is the total number of predicted values in a distribution (random sample size). The index becomes [(n) × (p ÷ 100)] and represents the place for the pth percentile value within the population distribution. All the predicted values are arranged from the smallest time value (in a week) until the largest time value in the distribution. For example, 0.581503 is the upper limit (95th percentile of N(t) at the end of 2015 for 5,000 simulations) of predicted values inside the population. It means 95 percent of predicted values for the population are equal to or less than 0.581503. The lower limit of this predicted interval is the 5th percentile. The median here is the 50th percentile of the predicted values. In other words, the 50th percentile is a data that split the entire data distribution into two pieces, such that 50 percent of the values lie in one piece and the other values lie in the second piece (Knowelden 1952). Here, we have chosen the 5th and 95th percentile to get a 90 percent prediction interval for our numerical-simulated solutions as an alternative to the 90 percent prediction interval of statistical predictions given by Santonja et al. (2012). Although the 5th and 95th percentiles are selected in the present work, we can also use other percentiles depending on the size of the prediction interval.

Results and Discussion

In this section, the MMCFD solutions are presented for the systems (1) to (3). Since random variables are included as parameters in this model, random variables sampling can be estimated by MMCFD method through MC simulation. Good inference is obtained when all MMCFD solutions are proved to be closer to statistical results (Santonja et al. 2012) than the FD numerical results with difference recurrence of simulations (100; 500; 1,000; 3,000; and 5,000 simulations) in 2013, 2015, and 2030 (Table 2). Since different approaches (numerical and statistical) are employed and due to source of errors, the MMCFD numerical-simulated solutions of the systems (1) to (3) are expected to have a percentage of difference with statistical and FD results spatially at the end of 2013, 2015, and 2030 for the subpopulations N(t), S(t), and O(t) as clearly supported by Tables 2 and 3 and Figures 2 to 4, respectively.

Table 3.

Absolute Growth Rates and Yearly Average Growth Rates of the Subpopulations.

Subpopulation	2000 Initial	Method	\|Growth\| (%)			\|Average Growth\| (%)
Subpopulation	2000 Initial	Method	2013	2015	2030	2013	2015	2030
N(t)	.704	From Santonja et al. (2012)	38.2	39.7	48.0	2.9	2.6	1.6
		FD	18.5	21.0	36.7	1.4	1.4	1.2
		MMCFD
		5,000 Simulation	19.9	22.9	42.6	1.5	1.5	1.4
S(t)	.250	From Santonja et al. (2012)	58.8	60.2	65.3	4.5	4.0	2.2
		FD	30.4	34.3	55.2	2.3	2.3	1.8
		MMCFD
		5,000 Simulation	32.1	36.6	61.6	2.5	2.4	2.1
O(t)	.046	From Santonja et al. (2012)	265.2	280.9	380.4	20.4	18.7	12.7
		FD	117.7	135.5	261.1	9.1	9.0	8.7
		MMCFD
		5,000 Simulation	126. 6	147.4	310.3	9.7	9.8	10.3

Figure 2.

Figure 3.

Figure 4.

In Table 2, solutions of N(t), S(t), and O(t) are tabulated in comparison with results from Santonja et al. (2012), FD method, and MMCFD method. The corresponding measure of differences between FD and MMCFD methods with respect to statistical predictions (Santonja et al. 2012) are also listed. In general, normal weight subpopulation, N(t), of the Valencia community is expected to reduce as the year increases from 2000 (see Table 1; N₀ = 0.704, S₀ = 0.250, O₀ = 0.046) to 2030, while the subpopulations of overweight S(t) and obese O(t) are growing over the same period. It is observed that the measure of difference of the MMCFD solutions consistently smaller than the measure of difference of the FD method in all years of 2013, 2015, and 2030 under various repeated simulations (100; 500; 1,000; 3,000; and 5,000 simulations). Although the current results depend on 5,000 simulations, different numbers of repeated simulations are tested in Table 2 to find the best number of simulation that produces the closest results to the predicted values (Santonja et al. 2012). In general, the statistical behavior is random with MMCFD method, and it is noted that MMCFD results are extremely convergent with regard to different number of simulations. It is observed that the MMCFD solutions for both normal weight and overweight subpopulations is the most convergent to the predicted values (Santonja et al. 2012) in 500 simulations in all years of 2013, 2015, and 2030. However, the MMCFD results of obese subpopulation in these years are the closest to the statistical values (Santonja et al. 2012) by using only 100 simulations.

On the other hand, the absolute growth rate percentage at the end of 2013, 2015, and 2030 as well as the absolute yearly average growth rate percentage as compared to the initial conditions of each subpopulation in the last week of 2000 (Santonja et al. 2012) are presented in Table 3. It is found that the normal weight N(t) and overweight S(t) subpopulations are decreasing at average rates of 1.6 percent and 2.2 percent per year by statistical prediction (Rohaninejad and Zarghami 2012), at 1.2 percent and 1.8 percent by FD method, and at 1.4 percent and 2.1 percent by MMCFD method, respectively, from the last week of 2000 until the end of 2030. Oppositely, the growth rate of obese subpopulation O(t) is increasing over the 30 years. The average yearly increasing rate for the obese subpopulation is estimated at 12.7 percent, 8.7 percent, and 10.3 percent from Santonja et al. (2012), FD, and MMCFD solutions, respectively. The statistical predictions (Santonja et al. 2012) produce the highest growth rates for all subpopulations and years, the FD method gives the minimum rates, while the MMCFD method provides the intermediate rates between the two different approaches of statistical and numerical methods.

According to Figure 2, normal weight curves for both FD and MMCFD methods are gradually declining. The decrement in N(t) is slightly higher in 2015 than in 2013. It is also noted that the normal weight subpopulations continue to fall until it approaches the predicted value (Santonja et al. 2012) at the end of 2030. This convinces the expectation that N(t) will decrease in future. In general, the normal weight curves for MMCFD are decreasing significantly than the normal weight curves of FD from 2001 to 2030. It is clear that MMCFD curves of normal weight subpopulation converge to the predicted values of Santonja et al. (2012) faster than the FD curves for all 30 years especially at the end of 2030.

Based on Figures 3 and 4, the overweight and obese curves are rising progressively in 2013. Also, the progress is more evidenced in 2015 as compared to 2013. Moreover, S(t) and O(t) are approaching the predicted values (Santonja et al. 2012) closely at the end of 2030. In other words, S(t) and O(t) are expected to increase with MMCFD overweight, and obese curves are consistently higher and closer to the predicted solutions (Santonja et al. 2012) than the FD curves in all years. Even though all curves of normal weight, overweight, and obese starting at the same initial points in the beginning of 2001, however only MMCFD solutions for all N(t), S(t), and O(t) are closer to the given predicted points (Santonja et al. 2012) at the end of 2013, 2015, and 2030 as compared to the FD solutions. It can be concluded that random error from simulation causes difference between FD and MMCFD results and the error increases with more iterations in the future as showed in Figures 2 to 4.

From Table 4, the standard deviation of mean, σ_n, in equations (17) to (20) of MMCFD for all subpopulations until 2030 is sufficiently small. Therefore, this indicator signifies an acceptable error (Jacoboni and Lugli 1989) for present work. From 2001 to 2030 at 5,000 simulations, standard deviation of error for normal weight subpopulation is not more than 0.006730, while the standard deviation of error for overweight individual drops to 0.005075. At the same time, the standard deviation of error for obese falls to 0.003191. The important point to be addressed here is the MMCFD standard deviation of mean, σ_n, has a very low percentage of error as evidenced in Table 4. Finally, σ_n increases with years as a result of statistical simulation errors.

Table 4.

Standard Deviation of Mean Monte Carlo Finite Difference Using 5,000 Simulations for Each Subpopulation.

Subpopulation	Standard Deviation of Mean, σ_n
Subpopulation	2013	2015	2030
N(t)	.002987	.003470	.006730
S(t)	.002789	.003171	.005075
O(t)	.001008	.001204	.003191

The predicted intervals of MMCFD solutions are approaching the range of intervals by Santonja et al. (2012) with a percentage of error (see Table 4 in Santonja et al. 2012 and Table 5). In general, many factors can cause errors in the present study. Firstly, modeling error might happen in the systems (1) to (3) as given in Santonja et al. (2012), where this error comes from prediction and actual data. Secondly, the sampling and rounding errors come from simulation method to estimate model parameters that are later being transported to predictions. The third source of the error is due to different simulation techniques to estimate model parameters. In the present work, the parameters are simulated by using the classical MC process, while in Santonja et al. (2012), the parameters are simulated 5,000 times by using LHS technique that is another type of MC sampling. Finally, the error comes from the type of distribution (uniform distribution) for parameters and solutions of the model (Santonja et al. 2012). The numerical method is subjected to two common numerical errors, the round off error which losing of precision due to computer rounding of decimal quantities and data propagation error. Therefore, the error is not restricted and extends in all interval areas of the distribution.

Table 5.

Prediction Intervals [5th Percentile, 95th Percentile] of Mean Monte Carlo Finite Difference Solutions for N(t), S(t), O(t) Using 5,000 Simulations at the End of 2013, 2015, and 2030.

Subpopulation	2013	2015	2030
N(t)	[.5297, .5970]	[.5030, .5815]	[.3226, .4806]
S(t)	[.2994, .3625]	[.3063, .3781]	[.3475, .4635]
O(t)	[.0931, .1162]	[.1004, .1283]	[.1542, .2291]

Because of random sampling and rounding errors in the particular model, the solutions can be presented by predicted interval of N(t), S(t), and O(t). The predicted interval consists of percentiles of final solutions for each variable N(t), S(t), and O(t) running at 5,000 evaluations, empirically at 5th and 95th percentiles. These percentiles are used to account 90 percent prediction interval that implies 90 percent of the predicted values inside this interval. Consider predicted intervals of distributions for N(t), S(t), and O(t) contain a lower limit (5th percentile) and an upper limit (95th percentile). Based on MMCFD solutions with 5,000 simulations, the predicted values inside the population is variously distributed for different years (see Table 5 and Figures 5 –7).

Figure 5.

Figure 6.

Figure 7.

Comparing between the present MMCFD results and the statistical predicted intervals, Santonja et al. (2012) presented 90 percent prediction intervals for N(t), S(t), and O(t) as listed in Table 4 (Santonja et al. 2012:37), which are taken from the empirical 5th and 95th percentiles. The current MMCFD-predicted interval consists of percentiles of N(t), S(t), and O(t) empirically at 5th and 95th percentiles too. Both predicted intervals of solutions to this model are obtained by using 5,000 simulations of model parameters. It is clear that all MMCFD solutions in Tables 2 and 6 with 5,000 simulations are inside the predicted intervals of percentile as presented in Table 5. In addition, the 50th percentile as the median for MMCFD solution distributions is close to the mean of these solutions as shown in all Figures 5 to 7 from 2001 to the end of 2013, 2015, and 2030, respectively.

Table 6.

The 50th Percentile (Median) of Mean Monte Carlo Finite Difference Results Using 5,000 Simulations With the Mean.

Subpopulation	Percentile	2013	2015	2030
N(t)	50th, Median	.564596	.543749	.405707
N(t)	Mean	.564039	.543108	.404293
S(t)	50th, Median	.329789	.340697	.403530
S(t)	Mean	.330336	.341476	.404109
O(t)	50th, Median	.103911	.113443	.186972
O(t)	Mean	.104218	.113818	.188732

Box plots analyzing the results of N(t), S(t), and O(t) from 2001 to 2030 are provided graphically in Figure 8. They show the central point (median), the extreme values (maximum and minimum values), and the outliers which are the points outside the boundaries. The box plots present the pattern of the MMCFD solutions of N(t), S(t), and O(t) throughout the 30 years. Using 5,000 simulations, the box plot of N(t) lays above in 2013 with one outlier point at bottom, declines significantly from 2015-2030. It means the subpopulation values of N(t) take smaller values in the end of 2030 as compared to 2013 and 2015 (see median in Table 6). Inversely, the box plots of S(t) and O(t) are located at the bottom in 2013, then they rise slightly in 2015 with major hikes in 2030. There is only one outlier value above the box for S(t) in 2030, but more outliers are at the top part of the box for O(t) in 2030. It means the subpopulation values of S(t) and O(t) have greater values at the end of 2030 than in 2013 and 2015 (see median in Table 6). The statistical behaviors shown in Figures 2 to 4 and in Figure 8 confirmed that N(t) will decrease while S(t) and O(t) will increase in the future.

Figure 8.

Box plots along 30 years (2001–2030) of mean Monte Carlo finite difference solutions With 5,000 simulations for (a) N(t), (b) S(t), (c) O(t).

Conclusion

In this study, the MMCFD method is proposed for the first time to solve an epidemiology model explicating the effects of public health campaign on body weight loss in Spain. The MMCFD method is suggested to generate new numerical-simulated solutions for this model, so that random variables sampling can be estimated by MC method. It is found that the MMCFD method produced numerical-simulated solutions of the model with better approximation toward previous statistical predictions as compared to the FD method in all time intervals considered. Through this work, it is shown that the MMCFD method is beneficial to generate random coefficient sampling for the nonlinear epidemic system. Besides, the MMCFD method is also promising to create agreement and estimation balance between a statistical method and a classical numerical method.

Footnotes

Acknowledgments

The authors thank the reviewers for their valuable suggestions and comments for the improvement of this article.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Abdulle

Blumenthal

. 2013. “Stabilized Multilevel Monte Carlo Method for Stiff Stochastic Differential Equations.” Journal of Computational Physics 251:445–60.

Buezas

F. S.

Rosales

M. B.

Sampaio

. 2013. “Propagation of Uncertainties and Multimodality in the Impact Problem of Two Elastic Bodies.” International Journal of Mechanical Sciences 75:145–55.

Christakis

N. A.

Fowler

J. H.

. 2007. “The Spread of Obesity in a Large Social Network Over 32 Years.” New England Journal of Medicine 357:370–79.

Dereich

Heidenreich

. 2011. “A Multilevel Monte Carlo Algorithm for Lévy-driven Stochastic Differential Equations.” Stochastic Processes and their Applications 121:1565–87.

De Souza

M. V. C.

Colaço

M. J.

Leiroz

A. J. K.

. 2014. “Application of the Generalized Polynomial Chaos Expansion to the Simulation of An Internal Combustion Engine with Uncertainties.” Fuel 134:358–67.

Jacoboni

Lugli

. 1989. The Monte Carlo Method for Semiconductor Device Simulation. 1st Ed. New York, NY: Springer.

Karris

S. T.

2007. Numerical Analysis Using MATLAB and Excel. 2nd Ed. Fremont, CA: Orchard. ISBN-13:978-1-934404-04-1.

Mohyud-Din

S. T.

Ali

. 2016. “On Biological Population Model of Fractional Order. International Journal of Biomathematics 9(5):1–13 (ID 1650070). doi:10.1142/s1793524516500704.

Knowelden

1952. “Introduction to Statistics.” British Medical Journal 2:270.

10.

Kovtanyuk

A. E.

Botkin

D. N.

Hoffmann

. 2012. “Numerical Simulations of a Coupled Radiative–Conductive Heat Transfer Model Using a Modified Monte Carlo Method.” International Journal of Heat and Mass Transfer 55:649–54.

11.

Leander

Lundh

Jirstrand

. 2014. “Stochastic Differential Equations as a Tool to Regularize the Parameter Estimation Problem for Continuous Time Dynamical Systems Given Discrete Time Measurements.” Mathematical Biosciences 251:54–62.

12.

Licea

J. A.

Villafuerte

Chen-Charpentier

B. M.

. 2013. “Analytic and Numerical Solutions of a Riccati Differential Equation with Random Coefficients.” Journal of Computational and Applied Mathematics 239:208–19.

13.

Livingston

E. H.

Sebastian

J. L.

Huerta

Yip

Heber

. 2001. “Biexponential Model for Predicting Weight Loss after Gastric Surgery for Obesity.” Journal of Surgical Research 101:216–24.

14.

Liu

Chen

. 2014. “Maximum Norm a Posteriori Error Estimates for a Singularly Perturbed Differential Difference Equation with Small Delay.” Applied Mathematics and Computation 227:801–10.

15.

Moghaddam

B. P.

Mostaghim

Z. S.

. 2013. “A Numerical Method Based on Finite Difference for Solving Fractional Delay Differential Equations. Journal of Taibah University for Science 7:120–27.

16.

Mohyud-Din

S. T.

Noor

M. A.

Noor

K. I.

. 2009. “Some Relatively New Techniques for Nonlinear Problems.” Mathematical Problems in Engineering 2009:1–25 (ID 234849). doi:10.1155/2009/234849.

17.

Mohyud-Din

S. T.

Yildirim

Gulkanat

. 2012. “Approximate Analysis of Population Dynamics with Density-dependent Migrations and the Allee Effects. International Journal of Numerical Methods for Heat & Fluid Flow 22:243–50.

18.

Navarro-Barrientos

Rivera

D. E.

Collins

L. M.

. 2011. “A Dynamical Model for Describing Behavioural Interventions for Weight Loss and Body Composition Change.” Mathematical and Computer Modelling of Dynamical Systems 17:183–203.

19.

Rao

R. N.

Chakravarthy

P. P.

. 2013. “A Finite Difference Method for Singularly Perturbed Differential-difference Equations Arising from a Model of Neuronal Variability. Journal of Taibah University for Science 7:128–36.

20.

Rao

R. N.

Chakravarthy

P. P.

. 2014. “An Exponentially Fitted Tridiagonal Finite Difference Method for Singularly Perturbed Differential-difference Equations with Small Shift. Ain Shams Engineering Journal 5:1351–60.

21.

Rohaninejad

Zarghami

. 2012. “Combining Monte Carlo and Finite Difference Methods for Effective Simulation of Dam Behavior.” Advances in Engineering Software 45:197–202.

22.

Rubinstein

R. Y.

Kroese

D. P.

. 2007. Simulation and the Monte Carlo Method. 2nd Ed. Hoboken, NJ: John Wiley and Sons.

23.

Santonja

F.-J.

Morales

Villanueva

R.-J.

Cortés

J.-C.

. 2012. “Analysing the Effect of Public Health Campaigns on Reducing Excess Weight: A Modeling Approach for the Spanish Autonomous Region of the Community of Valencia.” Evaluation and Program Planning. 35:34–39.

24.

Snyder

L. B.

2007. “Health Communication Campaigns and their Impact on Behavior.” Journal of Nutrition Education and Behavior 39:S32–40.

25.

Spigler

Zerbetto

. 2013. “Random Perturbations Effects on the Stability of Tension Leg Platforms in Offshore Engineering and of Other Large Structures.” Applied Mathematical Modelling. 37:2881–99.

26.

Ding

. 2013. “Numerical Dynamics of a Nonstandard Finite Difference Method for a Class of Delay Differential Equations. Journal of Mathematical Analysis and Applications 400:25–34.

27.

Thomas

D. M.

Weedermann

Fuemmeler

B. F.

Martin

C. K.

Dhurandhar

N. V.

Bredlau

Heymsfield

S. B.

Ravussin

Bouchard

. 2014. “Dynamic Model Predicting Overweight, Obesity, and Extreme Obesity Prevalence Trends.” Obesity 22:590–97.

28.

Ul Hassan

Q. M.

Mohyud-Din

S. T.

. 2016. “Investigating Biological Population Model Using Exp-function Method. International Journal of Biomathematics 9(2):1–15 (ID 1650026). doi:10.1142/s1793524516500261.

29.

Zhao

2013. “Compact Finite Difference Methods for High Order Integro-differential Equations. Applied Mathematics and Computation 221:66–78.