# Bayesian Analysis for a Fractional Population Growth Model

Authors: Francisco J. Ariza-Hernandez; Jorge Sanchez-Ortiz; Martin P. Arciga-Alejandre; Luis X. Vivas-Cruz
Journal: Journal of Applied Mathematics
 (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.


---

## Body

## 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 Section2. 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:(1)DαXt=aXt+ft,a>0,α∈0,1;(2)Xt0=X0,whereX(t) is the size of the population at the time t, X0 is the initial size, a is the rate of change, and f(t) is a harvesting function. Here Dα is the ordinary derivative for α=1 and the Caputo fractional derivative for α∈(0,1), which is defined by (3)DαXt=1Γ1-α∫0tX′τt-ταdτ,whereΓ is the Gamma function.(2)
An observation equation is as follows:(4)yi=hXti+εi,i=1,…,n,whereyi correspond to the ith observed value under uncertainty from a solution of (1) at the discrete time ti∈[0,T],i=1,2,…,n; h is the observation function; and εi 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 εi~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 (5)X^s=∫0∞e-stXtdt,Xt=12πi∫s0-i∞s0+i∞estX^sds,respectively. Applying the direct transform to (1), we obtain(6)X^s=1sα-asα-1X0+f^s.Then, using the following property (7)tβ-1Eα,βatα^s=sα-βsα-a,where Eα,β is the Mittag-Leffler function,(8)Eα,βz=∑k=0∞zkΓαk+β,α>0,β>0,z∈C.Finally, applying the inverse Laplace transform to (6) we arrive to(9)Xt=X0Eα,1atα+∫0tt-ηα-1Eα,αat-ηαfηdη.Note that ifα=1 in the above equation we get the well known solution for (1), which is given by(10)Xt=X0eat+∫0teat-ηfηdη.


## 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 θ′=(a,α,σ2,X0). Here, we propose the following prior distributions:(11)a~Gra,λa,(12)α~U0,1,(13)τ~Grτ,λτ,(14)X0~ln⁡Nμ0,σ02,where G(r,λ) denotes the Gamma distribution with shape parameter r and rate parameter λ, U is the continuous Uniform distribution on interval (0,1), τ=1/σ2 is called the precision parameter, and ln⁡N(μ0,σ02) is a log-Normal distribution with parameters μ0 and σ02. 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), a, X0, 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 (15)pθ∣hyperparameters∝pa∣ra,λapαpτ∣rτ,λτpX0∣μ0,σ02,where pa∣ra,λa, p(α), pτ∣rτ,λτ, and pX0∣μ0,σ02 are defined in (11)–(14), respectively. Let y′=(y1,y2,…,yn) denote i.i.d. observed data at times (t1,t2,…,tn) from the model defined by (1) and (4), and then the likelihood function is given by(16)Ly∣θ=τn/22π-n/2exp-τ2∑i=1nyi-hXti2.Thus, by Bayes’ theorem the posterior distribution of the parameters of interest is given by(17)pθ∣y=Ly∣θpθ∫ΘLy∣θpθdθ,where Θ denotes the parameter space of θ. It is clear that pθ∣y∝Ly∣θp(θ), and then the posterior distribution can be expressed by(18)pθ∣y∝τn/2exp-τ2∑i=1nyi-X0Eα,1atiα-∫0titi-ηα-1Eα,αati-ηαfηdη2ara-1τrτ-1X0-1exp-ln⁡X0-μ022σ02-aλa-τλτ.Using a loss quadratic function, the Bayesian point estimation is the posterior mean ofθ, which is given by (19)θ^B=Eθ∣y.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.Letf(θ) be univariate probability density function from which we want to obtain samples; for an initial value θ0, f(θ0)>0, the new value, θ1, will be found by a three-step procedure:(a)
Sample a real valuey~U(0,f(θ0)) and we define a horizontalslice, S={θ:y<f(θ)}.(b)
Find an interval,(l,r), 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θi,i=1,…,p, 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,M1, is given by the dynamical system (20)ddtXt=aXt,a>0;(21)Xt0=X0,and its corresponding observation equation, for n observations,(22)yi1=X0eati+εi,i=1,…,n.The fractional population growth model,M2, is given by(23)DαXt=aXt,a>0,α∈0,1;(24)Xt0=X0,(25)yi2=X0Eα,1atiα+εi,i=1,…,n,where Eα,1(·) is defined in (8) and εi~N(0,σ2). Note that for simplicity, in both models, the harvesting function f(t) is assumed to be zero and the function h is taken as the identity function.So, for modelMk, k=1,2, the posterior distribution takes the form(26)pθk∣y,Mk∝Ly∣θk,Mkpθk∣Mk,where Ly∣θk,Mk is the likelihood function of the data y′=y1,…,yn under model Mk and pθk∣Mk is the prior distribution of θk. By substituting these expressions for each model, we have (27)pθ1∣y,M1∝τn/2+rτ-1ara-1X0-1exp⁡-τ2∑i=1nyi-X0eati2-ln⁡X0-μ022σ02-λaa-λττ,pθ2∣y,M2∝τn/2+rτ-1ara-1X0-1exp-τ2∑i=1nyi-X0Eαatiα2-ln⁡X0-μ022σ02-λaa-λττ.To compare these models, we need to calculate the marginal likelihoodmy∣Mk for k=1,2, given by (28)my∣Mk=∫ΘkLy∣θk,Mkpθk∣Mkdθk,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 (29)my∣Mk=Ly∣θk,Mkpθk∣Mkpθk∣y,Mk.In general,my∣Mk 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 my∣Mk using (29) is not an easy task because we have many MCMC draws falling in the tails of pθk∣y,Mk and this leads to a very unreliable estimate of pθk∣y,Mk. Methods to estimate my∣Mk via MCMC sampling are discussed in [30–32]. We propose the use of Gelfand and Dey’s estimator [33] to obtain an approximation of my∣Mk, which is given by (30)m^GD=1G∑g=1GπθkgLy∣θkg,Mkpθkg∣Mk-1,where {θk(1),…,θk(G)} are draws from the posterior density pθk∣y,Mk 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 t densities with posterior sampling mean and covariance. Then, we compare the models M1 and M2 via Bayes factor:(31)BF1,2=exp⁡ln⁡my∣M1-ln⁡my∣M2,where a value of BF1,2 greater than one indicates a preference for model M1 with respect to the model M2.


## 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 θ′=(a,α,σ2,X0). Here, we propose the following prior distributions:(11)a~Gra,λa,(12)α~U0,1,(13)τ~Grτ,λτ,(14)X0~ln⁡Nμ0,σ02,where G(r,λ) denotes the Gamma distribution with shape parameter r and rate parameter λ, U is the continuous Uniform distribution on interval (0,1), τ=1/σ2 is called the precision parameter, and ln⁡N(μ0,σ02) is a log-Normal distribution with parameters μ0 and σ02. 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), a, X0, 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 (15)pθ∣hyperparameters∝pa∣ra,λapαpτ∣rτ,λτpX0∣μ0,σ02,where pa∣ra,λa, p(α), pτ∣rτ,λτ, and pX0∣μ0,σ02 are defined in (11)–(14), respectively. Let y′=(y1,y2,…,yn) denote i.i.d. observed data at times (t1,t2,…,tn) from the model defined by (1) and (4), and then the likelihood function is given by(16)Ly∣θ=τn/22π-n/2exp-τ2∑i=1nyi-hXti2.Thus, by Bayes’ theorem the posterior distribution of the parameters of interest is given by(17)pθ∣y=Ly∣θpθ∫ΘLy∣θpθdθ,where Θ denotes the parameter space of θ. It is clear that pθ∣y∝Ly∣θp(θ), and then the posterior distribution can be expressed by(18)pθ∣y∝τn/2exp-τ2∑i=1nyi-X0Eα,1atiα-∫0titi-ηα-1Eα,αati-ηαfηdη2ara-1τrτ-1X0-1exp-ln⁡X0-μ022σ02-aλa-τλτ.Using a loss quadratic function, the Bayesian point estimation is the posterior mean ofθ, which is given by (19)θ^B=Eθ∣y.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.Letf(θ) be univariate probability density function from which we want to obtain samples; for an initial value θ0, f(θ0)>0, the new value, θ1, will be found by a three-step procedure:(a)
Sample a real valuey~U(0,f(θ0)) and we define a horizontalslice, S={θ:y<f(θ)}.(b)
Find an interval,(l,r), 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θi,i=1,…,p, 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,M1, is given by the dynamical system (20)ddtXt=aXt,a>0;(21)Xt0=X0,and its corresponding observation equation, for n observations,(22)yi1=X0eati+εi,i=1,…,n.The fractional population growth model,M2, is given by(23)DαXt=aXt,a>0,α∈0,1;(24)Xt0=X0,(25)yi2=X0Eα,1atiα+εi,i=1,…,n,where Eα,1(·) is defined in (8) and εi~N(0,σ2). Note that for simplicity, in both models, the harvesting function f(t) is assumed to be zero and the function h is taken as the identity function.So, for modelMk, k=1,2, the posterior distribution takes the form(26)pθk∣y,Mk∝Ly∣θk,Mkpθk∣Mk,where Ly∣θk,Mk is the likelihood function of the data y′=y1,…,yn under model Mk and pθk∣Mk is the prior distribution of θk. By substituting these expressions for each model, we have (27)pθ1∣y,M1∝τn/2+rτ-1ara-1X0-1exp⁡-τ2∑i=1nyi-X0eati2-ln⁡X0-μ022σ02-λaa-λττ,pθ2∣y,M2∝τn/2+rτ-1ara-1X0-1exp-τ2∑i=1nyi-X0Eαatiα2-ln⁡X0-μ022σ02-λaa-λττ.To compare these models, we need to calculate the marginal likelihoodmy∣Mk for k=1,2, given by (28)my∣Mk=∫ΘkLy∣θk,Mkpθk∣Mkdθk,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 (29)my∣Mk=Ly∣θk,Mkpθk∣Mkpθk∣y,Mk.In general,my∣Mk 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 my∣Mk using (29) is not an easy task because we have many MCMC draws falling in the tails of pθk∣y,Mk and this leads to a very unreliable estimate of pθk∣y,Mk. Methods to estimate my∣Mk via MCMC sampling are discussed in [30–32]. We propose the use of Gelfand and Dey’s estimator [33] to obtain an approximation of my∣Mk, which is given by (30)m^GD=1G∑g=1GπθkgLy∣θkg,Mkpθkg∣Mk-1,where {θk(1),…,θk(G)} are draws from the posterior density pθk∣y,Mk 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 t densities with posterior sampling mean and covariance. Then, we compare the models M1 and M2 via Bayes factor:(31)BF1,2=exp⁡ln⁡my∣M1-ln⁡my∣M2,where a value of BF1,2 greater than one indicates a preference for model M1 with respect to the model M2.


## 4. Simulation Study

We carried out a simulation study to assess the performance of the Bayesian estimators of the parameters of interest (θk, k=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 simulatedM=1000 synthetic data sets form a fractional population growth model using the true values of the parameters α=0.5, a=0.2, σ2=30, and X0=20 for two sample sizes, n=10 and n=30, and two cases: X0 known and X0 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 modelM1 and for model M2.
(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 (32)a~G0.001,0.001,α~U0,1,τ~G0.001,0.001,note that for the parametersa and τ we use diffuse Gamma distributions, for α we propose Uniform prior distribution in (0,1), and, for the case of X0 unknown, a diffuse log-Normal distribution is used, X0~ln⁡N(0,1×103). These distributions are calledweakly 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 estimated95% 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 pθk∣y,Mk, k=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 ofM2versusM1, denoted by BF^2,1 using (31). If BF^2,1≥1, then we prefer the model M2 with respect to the model M1.(3)
From theseM 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 (n=10 and n=30) and two cases: X0 known and X0 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 M2 was more accurate compared to the model M1, in terms of the estimated bias and MSE. Also, the estimates coverage was near to 100%, in the most cases considered, for the model M2. From our simulation, for the case X0 known, the percentage of the estimated Bayes factors in favour of M2 was 98.4% with n=10 observations and 100% for n=30. For the case X0 unknown, the percentage of the estimated Bayes factors in favour of M2 was 94.0% with n=10 observations and 100% for n=30.Table 1
Estimated Bias, MSE, and coverage, withX0 known.
Parameter M 1 M 2 Bias MSE Cov (%) Bias MSE Cov (%) n = 10  a −0.105 0.011 0.0 0.008 0.002 97.4 α — — — 0.013 0.018 98.3 σ 2 25.22 1218.77 82.7 10.51 480.56 97.0  n = 30  a −0.135 0.018 0.0 −0.001 0.0004 92.3 α — — — 0.008 0.0019 92.4 σ 2 77.36 6457.78 0.0 1.091 71.55 94.0Table 2
Estimated Bias, MSE, and coverage, withX0 unknown.
Parameter M 1 M 2 Bias MSE Cov (%) Bias MSE Cov (%) n = 10  a −0.100 0.010 0.0 0.028 0.003 100 α — — — −0.004 0.005 100 σ 2 10.68 606.67 93.0 10.24 610.42 95.0 X 0 5.82 40.76 60.2 0.68 9.61 100  n = 30  a −0.151 0.023 0.0 −0.008 0.007 90.2 α — — — 0.092 0.029 90.1 σ 2 5.97 131.63 90.1 1.60 79.42 95.3 X 0 9.44 90.63 0.0 1.44 14.69 92.0


## 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 Figure2, 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 α, a, and σ2.The prior distributions used were forα~U(0,1), a~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.Figure 1
Trace and estimated posterior densities ofa, α, and σ2.Figure 2
Mobile web (in percentage) from December 2008 to August 2014.Posterior means are used as point estimates of the parameters of interest, forM1:a^=0.06 and σ^2=4.31 and for M2:α^=0.35, a^=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 M2.


## 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,α=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.



---

*Source: 9654506-2017-01-23.xml*