Bayesian Analysis for a Fractional Population Growth Model

We implement the Bayesian statistical inversion theory to obtain the solution for an inverse problem of growth data, using a fractional population growth model. We estimate the parameters in the model and we make a comparison between this model and an exponential one, based on an approximation of Bayes factor. A simulation study is carried out to show the performance of the estimators and the Bayes factor. Finally, we present a real data example to illustrate the effectiveness of the method proposed here and the pertinence of using a fractional model.


Introduction
In recent years, models governed by fractional differential equations have been of great importance to describe anomalous dynamic, which is inherent in a wide range of phenomena (see [1][2][3]).This model under certain circumstances can offer a superior fit to experimental data [4].It is well known that population growth can be modeled by ordinary differential equations and the exponential growth model is the simplest and most widely used of all.However, it can fail to describe an anomalous growth; that is, the population may have a slower or faster growth than exponential growth model.In order to fill this gap, we use a fractional population growth model.
On the order hand, the measurement procedure of growth data is under uncertainty or error data collection process.The statistical approach toward inverse problems from a Bayesian point of view has been very important for the development of uncertainty quantification.Textbooks such as [5][6][7] provide an extensive literature on the Bayesian statistical aspects in inverse problems.Donnet et al. [8] developed a Bayesian inference for the parameters of stochastic differential equations deduced from the standard deterministic growth function by adding random effects to the growth dynamic in a population of individuals over time.Capistrán et al. [9] also performed a Bayesian analysis of inverse problems in ordinary differential equations.Calvetti et al. [10] mentioned several articles on the development of Bayesian inverse problems in a broad range of applications, such as engineering, geophysics, life sciences, and economy.Dashti and Stuart [11] discussed the Bayesian approach to inverse problems in differential equations from mathematical and computational perspective and describe Markov Chain Monte Carlo and sequential Monte Carlo methods and measure-preserving reversible stochastic differential equations on infinite dimensional space.
In this paper, we implement the Bayesian statistical inversion theory to obtain a solution for an inverse problem of growth data using a fractional population growth model, defined in Section 2. In Section 3, we estimate the parameters of the fractional model proposed and we make a comparison between two different models, based on Bayes factor.In Section 4, we make a simulation study to show the performance of the estimators and the Bayes factor.Finally, we present a real data example to illustrate the effectiveness of the method proposed here and the pertinence of using a fractional model.

The Model
In this work, the proposed mathematical model is defined by two parts.

Journal of Applied Mathematics
(1) A dynamical system is as follows: D   () =  () +  () ,  > 0,  ∈ (0, 1] ; (1) where () is the size of the population at the time ,  0 is the initial size,  is the rate of change, and () is a harvesting function.Here D  is the ordinary derivative for  = 1 and the Caputo fractional derivative for  ∈ (0, 1), which is defined by where Γ is the Gamma function.(2) An observation equation is as follows: where   correspond to the th observed value under uncertainty from a solution of (1) at the discrete time   ∈ [0, ],  = 1, 2, . . ., ; ℎ is the observation function; and   are measurement errors, which are considered as independent and identically distributed (i.i.d.) random variables from a normal distribution, with mean zero and constant variance  2 , denoted by   ∼ N(0,  2 ) (in this paper, the symbol "∼" is used for "is distributed as").
Now, we find an explicit solution for (1) using the Laplace transform and its inverse, which are defined by respectively.Applying the direct transform to (1), we obtain Then, using the following property where  , is the Mittag-Leffler function, Finally, applying the inverse Laplace transform to (6) we arrive to Note that if  = 1 in the above equation we get the well known solution for (1), which is given by

Bayesian Analysis
3.1.Bayesian Estimation.In recent years, the problem of parameter estimation in fractional models has been of great interest, and we use a Bayesian approach to estimate the parameters in a fractional population growth model.The Bayesian theory has proven to be efficient to solve fractional inverse problems (see [12,13]).In this setting, the parameters are regarded as random variables, and they have prior distributions that reflect the knowledge about their true values before observing the data.From the perspective of Bayesian statistical inversion theory, the solution to an inverse problem is the posterior distribution of the quantity of interest given that all information available has been incorporated in the model [5].
Note that, for the model defined in ( 1) and ( 4), the parameter of interest is   = (, ,  2 ,  0 ).Here, we propose the following prior distributions: ∼ U (0, 1) , where G(, ) denotes the Gamma distribution with shape parameter  and rate parameter , U is the continuous Uniform distribution on interval (0, 1),  = 1/ 2 is called the precision parameter, and ln N( 0 ,  2 0 ) is a log-Normal distribution with parameters  0 and  2 0 .The parameters involved in the prior distribution are called hyperparameters.The prior distributions were chosen according to our knowledge about the possible values of the parameters of interest.In models (1) and ( 4), ,  0 , and  are positive and  take values in the interval (0, 1].Assuming prior independence of the parameters, we can write the joint prior distribution as where ( |   ,   ), (), ( |   ,   ), and ( 0 |  0 ,  2 0 ) are defined in ( 11)-( 14), respectively.Let y  = ( 1 ,  2 , . . .,   ) denote i.i.d.observed data at times ( 1 ,  2 , . . .,   ) from the model defined by ( 1) and ( 4), and then the likelihood function is given by Thus, by Bayes' theorem the posterior distribution of the parameters of interest is given by where Θ denotes the parameter space of .It is clear that ( | y) ∝ (y | )(), and then the posterior distribution can be expressed by Using a loss quadratic function, the Bayesian point estimation is the posterior mean of , which is given by Note that the joint posterior distribution is analytically intractable and the marginal posterior distributions of the parameters are complicated; however, we can obtain samples from (18) using Markov Chain Monte Carlo (MCMC) techniques.The most common MCMC methods are the Gibbs sampling [14,15] and the Metropolis-Hastings [16,17].Currently, many of the MCMC algorithms have been already implemented in computer programs, such as WinBUGS and JAGS [18] and Stan and t-walk [19].All of these software packages provide programs for Bayesian modeling through posterior simulation given a specified model and data.The R [20] packages, such as R2WinBUGS [21], R2jags [22], rjags [23], and rstan [24], allow one to run WinBUGS, JAGS, and Stan from within R software, respectively.In this paper, we use JAGS within R to obtain samples from the marginal posterior distributions of interest.
JAGS (Just Another Gibbs Sampler) is a program for analysis of Bayesian models using MCMC methods, which was written in C++ by Martyn Plummer with three objectives: (a) to have a cross-platform engine for the BUGS language, (b) to be extensible, allowing users to write their own functions, distributions, and samplers, and (c) to be a platform for experimentation with ideas in Bayesian modeling (see [18]); moreover, it is a free package.To draw samples from the posterior distribution, JAGS use the Gibbs sampling and sometimes it is combined with or replaced by complex techniques such as Metropolis-Hastings algorithm, slice sampling [25], and Adaptive Rejection sampling [26].In our work, JAGS used slice sampling, which is briefly described as follows.
Let () be univariate probability density function from which we want to obtain samples; for an initial value  0 , ( 0 ) > 0, the new value,  1 , will be found by a three-step procedure: (a) Sample a real value  ∼ U(0, ( 0 )) and we define a horizontal slice,  = { :  < ()}.(b) Find an interval, (, ), around  0 that contains all, or much, of the slice.
(c) Sample the new value,  1 , from the part of the slice within this interval.
For multivariate distributions, we can apply the single variable slice sampling for each   ,  = 1, . . ., , in turn.An example of JAGS implementation is given in Appendix.

Bayesian Model Comparison.
In this section, our interest is to make a comparison between two different models.In one, the dynamic system changes in an ordinary way and in the other it changes anomalously, that is, with an entire derivative and with a fractional derivative, respectively.As soon as we obtain experimental data from a dynamical phenomenon, the uncertainty becomes present.Moreover, the mechanism that governs the evolution law is unknown at some level.So, the challenge here is to know which model best fits and best describes the data.
The Bayesian model selection in inverse problems framework has not been much investigated.In [8] it is suggested to validate the SDE approach via criteria based on the predictive posterior distribution.To compare a numerical and the theoretical posterior distribution, in [9] it was proposed to use Bayes factors, considering both of them as models for the data at hand.Other authors (see [27][28][29]) use Bayes factor to select the best model.We carried out a comparison of two models, defined below, from Bayesian point of view via Bayes factor, where the marginal likelihood is obtained using Gelfand and Dey's estimator.Now, let us define the models that we are going to compare.The ordinary population growth model, M 1 , is given by the dynamical system    () =  () ,  > 0; (20) and its corresponding observation equation, for  observations, The fractional population growth model, M 2 , is given by where  ,1 (⋅) is defined in (8) and   ∼ N(0,  2 ).Note that for simplicity, in both models, the harvesting function () is assumed to be zero and the function ℎ is taken as the identity function.
So, for model M  ,  = 1, 2, the posterior distribution takes the form where (y |   , M  ) is the likelihood function of the data y  = ( 1 , . . .,   ) under model M  and (  | M  ) is the prior distribution of   .By substituting these expressions for each model, we have To compare these models, we need to calculate the marginal likelihood (y | M  ) for  = 1, 2, given by and we choose the model which yields the largest marginal likelihood [30].Note that the equation above is the normalizing constant of the posterior density in (26), so it can be written as In general, (y | M  ) in (28) can not be represented in terms of elementary functions and obtaining a numerical approximation usually requires a huge computational effort.On the other hand, to calculate (y | M  ) using ( 29) is not an easy task because we have many MCMC draws falling in the tails of (  | y, M  ) and this leads to a very unreliable estimate of (  | y, M  ).Methods to estimate (y | M  ) via MCMC sampling are discussed in [30][31][32].We propose the use of Gelfand and Dey's estimator [33] to obtain an approximation of (y | M  ), which is given by where { (1)   , . . .,  ()   } are draws from the posterior density (  | y, M  ) obtained using MCMC method and (⋅) is an importance sampling density.The stability of (30) depends on the importance density choice, which is recommended to be of thinner tails compared to the product of the prior and the likelihood.Natural choices of (⋅) would be the multivariate normal or Student  densities with posterior sampling mean and covariance.Then, we compare the models M 1 and M 2 via Bayes factor: where a value of BF 1,2 greater than one indicates a preference for model M 1 with respect to the model M 2 .

Simulation Study
We carried out a simulation study to assess the performance of the Bayesian estimators of the parameters of interest (  ,  = 1,2) and the performance of the Bayes factor using Gelfand and Dey's estimator [32] for the marginal likelihood (29).For this task, we used R software [20] with R2jags package [22], under the following algorithm: (1) We simulated  = 1000 synthetic data sets form a fractional population growth model using the true values of the parameters  = 0.5,  = 0.2,  2 = 30, and  0 = 20 for two sample sizes,  = 10 and  = 30, and two cases:  0 known and  0 unknown.Each data set was simulated using the model defined in (25).
(2) For each simulated sample, we obtained Bayesian estimators of the parameters both for model M 1 and for model M 2 .
(a) We calculated the Bayes estimators as the sample mean from ( 18) using JAGS within R. We used two chains, each with 10000 iterations with a burn-in of 1000 iterations and a thinning rate of 9, so we kept a total of 2000 iterations to make inferences about the parameters of interest.The following prior distributions for the parameters in the models are assumed as  ∼ G (0.001, 0.001) ,  ∼ U (0, 1) ,  ∼ G (0.001, 0.001) , (32) note that for the parameters  and  we use diffuse Gamma distributions, for  we propose Uniform prior distribution in (0, 1), and, for the case of  0 unknown, a diffuse log-Normal distribution is used,  0 ∼ ln N(0, 1×10 3 ).These distributions are called weakly informative prior distributions, which contain some information about the parameters of interest but without affecting inferences by information external to the current data [34].(b) We estimated 95% credible intervals for each parameter in each model; these intervals are based on the 0.025 quantile and the 0.975 quantile of the corresponding posterior sample.(c) We calculate the marginal likelihood using Gelfand and Dey's estimator given in (30) using MCMC samples from the posterior density (  | y, M  ),  = 1, 2. The multivariate normal density with posterior sampling mean and covariance was used as importance sampling density.(d) We obtain an estimation of the Bayes factor of M 2 versus M 1 , denoted by BF 2,1 using (31).If BF 2,1 ≥ 1, then we prefer the model M 2 with respect to the model M 1 .
(3) From these  samples, we obtained measures of the bias and mean-squared error of the estimators, as well as the coverage of the corresponding credible intervals.
Our simulations were carried out for two sample sizes ( = 10 and  = 30) and two cases:  0 known and  0 unknown.Tables 1 and 2 show the estimate bias, MSE, and coverage for the cases considered.We can observe that the parameters estimation of the model M 2 was more accurate compared to the model M 1 , in terms of the estimated bias and MSE.Also, the estimates coverage was near to 100%, in the most cases considered, for the model M 2 .From our simulation, for the case  0 known, the percentage of the estimated Bayes factors in favour of M 2 was 98.4% with  = 10 observations and 100% for  = 30.For the case  0 unknown, the percentage of the estimated Bayes factors in favour of M 2 was 94.0% with  = 10 observations and 100% for  = 30.

Application
In this section, an example is given to illustrate the Bayesian estimation and selection model approach proposed in this paper.The data set used consists of 69 monthly measurements (in percentages) of the mobile web use in the world, from December 2008 to August 2014 (see Figure 2, dot points).We can find this data from the link http://stats.areppim.com/stats/statsmobiwebsubstxtime.htm.At a first glance of data, one can think of an exponential growth model; however, as we will see, this model is unable to fit the data in satisfactory way.A generalization of this model is proposed, where we add a new parameter for the order of the derivative.The model defined in (23) is used to fit this data, and the parameters of interest to be estimated, based in available data, are , , and  2 .
The prior distributions used were for  ∼ U(0, 1),  ∼ G(0.001, 0.001), and 1/ 2 ∼ G(0.001, 0.001).For the model proposed, we used two chains, each with 500000 iterations, and the first 250000 were discarded, taking a thinning rate of 50.Thus, 10000 posterior samples were used to obtain the summary statistics about the parameters of interest.Standard convergence diagnostics were carried out.To mention a few, the value of R for each parameter of interest was quite close to 1 for all models considered.Also, the Gelman-Rubin factor and Geweke diagnostics were calculated and showed evidence of convergence.Figure 1 shows the trace and estimated posterior distributions of the estimated parameters of interest.library(R2jags) library(coda) library(lattice) library(R2WinBUGS) library(rjags) Posterior means are used as point estimates of the parameters of interest, for M 1 : â = 0.06 and σ2 = 4.31 and for M 2 : α = 0.35, â = 0.33, and σ2 = 0.52. Figure 2 shows the fitted models considered.The Bayes factor was 64.25, giving very strong evidence in favour of the model M 2 .

Concluding Remarks
We estimated the parameters of a fractional population growth model using Bayesian framework.From the model proposed we can distinguish two different models, one with ordinary derivative,  = 1, and another with fractional derivative,  ∈ (0, 1).An estimation of the Bayes factor was used to compare these two growth models.To measure the performance of the parameters estimation and Bayes factor a simulation study was conducted.As expected, the fractional population growth model outperforms the ordinary population growth model as is shown in Tables 1 and 2. In order to estimate the Bayes factor, Gelfand and Dey's estimator was used to approximate the marginal likelihood.It is known that this estimator could be unstable, even in the situation when the importance density is chosen as a multivariate normal density with mean and covariance obtained from MCMC sampling; however, in this work, it was not the case.A subsequent monitoring was carried on to check the stability.Our results give numerical evidence of good performance of the Bayes factor.Also, the inferences based on posterior simulation are relatively easy to implement in JAGS within the R software.Finally, we have shown a real data example where the fractional growth model had a better fit with respect to the ordinary model.

( 5 )
Initial Values.Define the starting values for the MCMC runs (see Box 5).

( 7 )
Diagnostic.Convert the model output into an MCMC object in order to have access to several convergence diagnostics (seeBox 7).