Bayesian Analysis for a Fractional Population Growth Model

Ariza-Hernandez, Francisco J.; Sanchez-Ortiz, Jorge; Arciga-Alejandre, Martin P.; Vivas-Cruz, Luis X.

doi:https://doi.org/10.1155/2017/9654506

Journal of Applied Mathematics

On this page

Abstract Introduction Appendix References Copyright Related Articles

Research Article | Open Access

Volume 2017 | Article ID 9654506 | https://doi.org/10.1155/2017/9654506

Bayesian Analysis for a Fractional Population Growth Model

Francisco J. Ariza-Hernandez,¹Jorge Sanchez-Ortiz,¹Martin P. Arciga-Alejandre,¹and Luis X. Vivas-Cruz¹

Academic Editor: Juan J. Trujillo

Received26 Jul 2016

Revised04 Nov 2016

Accepted27 Nov 2016

Published23 Jan 2017

Abstract

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.

1. 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–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–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.

2. The Model

In this work, the proposed mathematical model is defined by two parts.(1)A dynamical system is as follows: where is the size of the population at the time , is the initial size, is the rate of change, and is a harvesting function. Here is the ordinary derivative for and the Caputo fractional derivative for , 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 ; 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 , denoted by (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 obtainThen, using the following property where is the Mittag-Leffler function,Finally, applying the inverse Laplace transform to (6) we arrive to

Note that if in the above equation we get the well known solution for (1), which is given by

3. 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 . Here, we propose the following prior distributions:where denotes the Gamma distribution with shape parameter and rate parameter , is the continuous Uniform distribution on interval , is called the precision parameter, and is a log-Normal distribution with parameters and . 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), , , and are positive and take values in the interval . Assuming prior independence of the parameters, we can write the joint prior distribution as where , , , and are defined in (11)–(14), respectively. Let denote i.i.d. observed data at times 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 bywhere denotes the parameter space of . It is clear that , 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 , , the new value, , will be found by a three-step procedure:(a)Sample a real value and we define a horizontal slice, .(b)Find an interval, , around that contains all, or much, of the slice.(c)Sample the new value, , from the part of the slice within this interval.

For multivariate distributions, we can apply the single variable slice sampling for each , in turn. An example of JAGS implementation is given in Appendix.

3.2. 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–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, , is given by the dynamical system and its corresponding observation equation, for observations,

The fractional population growth model, , is given bywhere is defined in (8) and . 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 , , the posterior distribution takes the formwhere is the likelihood function of the data under model and 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 for , 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, 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 using (29) is not an easy task because we have many MCMC draws falling in the tails of and this leads to a very unreliable estimate of . Methods to estimate via MCMC sampling are discussed in [30–32]. We propose the use of Gelfand and Dey’s estimator [33] to obtain an approximation of , which is given by where are draws from the posterior density 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 and via Bayes factor:where a value of greater than one indicates a preference for model with respect to the model .

4. Simulation Study

We carried out a simulation study to assess the performance of the Bayesian estimators of the parameters of interest (, ) 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 synthetic data sets form a fractional population growth model using the true values of the parameters , , , and for two sample sizes, and , and two cases: known and 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 and for model .(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 note that for the parameters and we use diffuse Gamma distributions, for we propose Uniform prior distribution in , and, for the case of unknown, a diffuse log-Normal distribution is used, . 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 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 , . 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 versus , denoted by using (31). If , then we prefer the model with respect to the model .(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 ( and ) and two cases: known and 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 was more accurate compared to the model , in terms of the estimated bias and MSE. Also, the estimates coverage was near to 100%, in the most cases considered, for the model . From our simulation, for the case known, the percentage of the estimated Bayes factors in favour of was 98.4% with observations and 100% for . For the case unknown, the percentage of the estimated Bayes factors in favour of was 94.0% with observations and for .

5. 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/stats_mobiwebsubstxtime.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 .

The prior distributions used were for , , and . 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 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.

Posterior means are used as point estimates of the parameters of interest, for and and for , , and . Figure 2 shows the fitted models considered. The Bayes factor was , giving very strong evidence in favour of the model .

6. 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, , and another with fractional derivative, . 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.

Appendix

JAGS Implementation

The proposed models were all implemented in JAGS using the R2jags package. R2jags was used to fit the models and to perform convergence diagnostics right within R. Here we use the data set of Section 5 to illustrate the implementation of model in JAGS.

(1) Packages. Load the required R packages (see Box 1).

(2) Data. Read the data from the working directory (see Box 2).

(3) The Model. Write the model in BUGS code and save it as “movil.model.jags” in the working directory (see Box 3).

(4) Parameters. Define the parameters of interest (see Box 4).

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

(6) Fit. Fit the model in JAGS (see Box 6).

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

Competing Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

C. A. Monje, Y. Chen, B. Vinagre, D. Xue, and V. Feliu, Fractional-Order Systems and Controls: Fundamentals and Applications, Advances in Industrial Control, Springer Science & Business Media, London, UK, 2010.
View at: Publisher Site | MathSciNet
L. Vázquez, J. J. Trujillo, and M. P. Velasco, “Fractional heat equation and the second law of thermodynamics,” Fractional Calculus and Applied Analysis, vol. 14, no. 3, pp. 334–342, 2011.
View at: Publisher Site | Google Scholar | MathSciNet
I. Area, H. Batarfi, J. Losada, J. J. Nieto, W. Shammakh, and Á. Torres, “On a fractional order Ebola epidemic model,” Advances in Difference Equations, vol. 2015, article 278, 2015.
View at: Publisher Site | Google Scholar
B. Jin and W. Rundell, “A tutorial on inverse problems for anomalous diffusion processes,” Inverse Problems. An International Journal on the Theory and Practice of Inverse Problems, Inverse Methods and Computerized Inversion of Data, vol. 31, no. 3, Article ID 035003, 2015.
View at: Publisher Site | Google Scholar | MathSciNet
J. Kaipio and E. Somersalo, Statistical and Computational Inverse Problems, vol. 160 of Applied Mathematical Sciences, Springer, New York, NY, USA, 2006.
View at: MathSciNet
H. T. Banks, S. Hu, and W. C. Thompson, Modeling and Inverse Problems in the Presence of Uncertainty, Monographs and Research Notes in Mathematics, CRC Press, Taylor & Francis, Boca Raton, Fla, USA, 2014.
View at: MathSciNet
A. Tarantola, Inverse Problem Theory and Methods for Model Parameter Estimation, SIAM, 2005.
View at: Publisher Site | MathSciNet
S. Donnet, J.-L. Foulley, and A. Samson, “Bayesian analysis of growth curves using mixed models defined by stochastic differential equations,” Biometrics, vol. 66, no. 3, pp. 733–741, 2010.
View at: Publisher Site | Google Scholar | MathSciNet
M. A. Capistrán, J. A. Christen, and S. Donnet, “Bayesian analysis of ODEs: solver optimal accuracy and Bayes factors,” SIAM/ASA Journal on Uncertainty Quantification, vol. 4, no. 1, pp. 829–849, 2016.
View at: Publisher Site | Google Scholar | MathSciNet
D. Calvetti, J. P. Kaipio, and E. Somersalo, “Inverse problems in the Bayesian framework,” Inverse Problems, vol. 30, no. 11, Article ID 110301, 2014.
View at: Publisher Site | Google Scholar | MathSciNet
M. Dashti and A. M. Stuart, “The Bayesian approach to inverse problems,” in Handbook of Uncertainty Quantification, Lecture Notes in Computer Science, Edited by R. Ghanem, D. Higdon and H. Owhadi, Springer, Berlin, Germany, 2016.
View at: Google Scholar
W. Fan, X. Jiang, and H. Qi, “Parameter estimation for the generalized fractional element network Zener model based on the Bayesian method,” Physica A. Statistical Mechanics and its Applications, vol. 427, pp. 40–49, 2015.
View at: Publisher Site | Google Scholar | MathSciNet
W. Fan, X. Jiang, and S. Chen, “Parameter estimation for the fractional fractal diffusion model based on its numerical solution,” Computers & Mathematics with Applications, vol. 71, no. 2, pp. 642–651, 2016.
View at: Publisher Site | Google Scholar | MathSciNet
A. E. Gelfand and A. F. M. Smith, “Sampling-based approaches to calculating marginal densities,” Journal of the American Statistical Association, vol. 85, no. 410, pp. 398–409, 1990.
View at: Publisher Site | Google Scholar | MathSciNet
G. Casella and E. I. George, “Explaining the Gibbs sampler,” The American Statistician, vol. 46, no. 3, pp. 167–174, 1992.
View at: Publisher Site | Google Scholar | MathSciNet
N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” The Journal of Chemical Physics, vol. 21, no. 6, pp. 1087–1092, 1953.
View at: Publisher Site | Google Scholar
S. Chib and E. Greenberg, “Understanding the metropolis-hastings algorithm,” The American Statistician, vol. 49, no. 4, pp. 327–335, 1995.
View at: Publisher Site | Google Scholar
M. Plummer, JAGS: A program for analysis of bayesian graphical models using Gibbs sampling, 2012.
J. Andrés Christen, C. Fox, D. Andrés Pérez-Ruiz, and M. Santana-Cibrian, “On optimal direction gibbs sampling,” https://arxiv.org/abs/1205.4062.
View at: Google Scholar
R Core Team, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, 2014.
S. Sturtz, U. Ligges, and A. Gelman, “R2WinBUGS: a package for running WinBUGS from R,” Journal of Statistical Software, vol. 12, no. 3, pp. 1–16, 2005.
View at: Google Scholar
Y.-S. Su and M. Yajima, R2jags: A Package for Running “JAGS’ from R, R package version 0.05-03, 2015.
M. Plummer, rjags: Bayesian Graphical Models using MCMC, R package version 3–15, 2015.
Stan Development Team, Stan Modeling Language User's Guide and Reference Manual, Version 2.6.1., 2015.
R. M. Neal, “Slice sampling,” The Annals of Statistics, vol. 31, no. 3, pp. 705–767, 2003.
View at: Publisher Site | Google Scholar | MathSciNet
W. R. Gilks and P. Wild, “Adaptive rejection sampling for gibbs sampling,” Applied Statistics, vol. 41, no. 2, pp. 337–348, 1992.
View at: Publisher Site | Google Scholar
Y. Murakami, “Bayesian parameter inference and model selection by population annealing in systems biology,” PLOS ONE, vol. 9, no. 8, Article ID e104057, 2014.
View at: Publisher Site | Google Scholar
X. Zhang, J. Cao, and R. J. Carroll, “On the selection of ordinary differential equation models with application to predator-prey dynamical models,” Biometrics, vol. 71, no. 1, pp. 131–138, 2015.
View at: Publisher Site | Google Scholar | MathSciNet
S. Hug, D. Schmidl, W. B. Li, M. B. Greiter, and F. J. Theis, “Uncertainty in biology,” in Bayesian Model Selection Methods and Their Application to Biological ODE Systems, pp. 243–268, Springer International, Geneva, Switzerland, 2016.
View at: Google Scholar
M.-H. Chen, Q.-M. Shao, and J. G. Ibrahim, Monte Carlo Methods in Bayesian Computation, Springer Science & Business Media, Berlin, Germany, 2012.
R. E. Kass and A. E. Raftery, “Bayes factors,” Journal of the American Statistical Association, vol. 90, no. 430, pp. 773–795, 1995.
View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
S. Chib, “Marginal likelihood from the Gibbs output,” Journal of the American Statistical Association, vol. 90, no. 432, pp. 1313–1321, 1995.
View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
A. E. Gelfand and D. K. Dey, “Bayesian model choice: asymptotics and exact calculations,” Journal of the Royal Statistical Society. Series B. Methodological, vol. 56, no. 3, pp. 501–514, 1994.
View at: Google Scholar | Zentralblatt MATH | MathSciNet
A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vethari, and D. B. Rubin, Bayesian Data Analysis, Taylor &Francis, 3rd edition, 2014.

Copyright

Copyright © 2017 Francisco J. Ariza-Hernandez et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

2129

Downloads

1172

Citations