Abstract
Extensive research has been devoted to the estimation of the parameters of frequently used distributions. However, little attention has been paid to estimation of parameters of Gamma/Gompertz distribution, which is often encountered in customer lifetime and mortality risks distribution literature. This distribution has three parameters. In this paper, we proposed an algorithm for estimating the parameters of Gamma/Gompertz distribution based on maximum likelihood estimation method. Iterated local search (ILS) is proposed to maximize likelihood function. Finally, the proposed approach is computationally tested using some numerical examples and results are analyzed.
1. Introduction
In probability and statistics, the Gompertz distribution is a continuous probability distribution. The Gompertz distribution was first introduced by Gompertz [1] to describe human mortality and establish actuarial tables. Since then, many investigators have used the Gompertz distribution or some related forms of it in a variety of studies. There are many forms of the Gompertz distribution in the literature. The Gompertz distribution is often applied to describe the distribution of adult lifespans by demographers [2, 3] and actuaries [4, 5]. Related fields of science such as biology [6] and gerontology [7] also considered the Gompertz distribution for the analysis of survival. More recently, computer scientists have also started to model the failure rates of computer codes by the Gompertz distribution [8]. In marketing science, it has been used as an individual-level model of customer lifetime [9]. Also, the Gompertz distribution is applied in reliability, life testing, epidemiological, and biomedical studies. Several such situations have been discussed by Ananda et al. [10], Walker and Adham [11], Jaheen [12], and many others. The probability density function of the Gompertz distribution is where is the scale parameter and is the shape parameter of the Gompertz distribution. In the actuarial and biological sciences and in demography, the Gompertz distribution is parameterized slightly and differently The Gompertz distribution is a flexible distribution that can be skewed to the right and to the left. The Gompertz density function can take on different shapes depending on the values of the shape parameter . The th moment of a Gompertz distributed random variable is where is the generalized integroexponential function [13]. When the shape parameter of a Gompertz distribution varies according to a Gamma distribution with shape parameter and scale parameter (), the distribution of is Gamma/Gompertz [9]. The Gamma/Gompertz distribution has been used as an aggregate-level model of customer lifetime and a model of mortality risks. Missov studied life expectancy resulting from a gamma-Gompertz force of mortality [14].
The probability density function of the Gamma/Gompertz distribution is where is the scale parameter and and are the shape parameters of the Gamma/Gompertz distribution (Figure 1). The cumulative distribution function of the Gamma/Gompertz distribution is When , this reduces to an exponential distribution with parameter : The expected value of random variable of a Gamma/Gompertz distribution depends on the values of , , and ; that is, we have for , and for , . Also, the median of this distribution is .
Successful application of Gamma/Gompertz distribution depends on having acceptable statistical estimates of the distribution parameters. One of the best methods of obtaining a point estimator of a parameter is the method of maximum likelihood. As the name implies, the estimator will be the value of the parameter that maximizes the likelihood function. While the method of maximum likelihood is an excellent technique, sometimes complications arise in its use. For example, it may not always be possible to use calculus methods directly to determine the maximum of likelihood function. This difficulty is obvious in estimating the parameters of the Gamma/Gompertz distribution.
In this paper, we use iterated local search (ILS) to estimate Gamma/Gompertz parameters. The paper is organized as follows. Section 2 describes the notion of maximum likelihood estimation. Section 3 explains the basics of iterated local search (ILS). In Section 4, we explain the steps of ILS algorithm to solve the problem. Some numerical examples and their computational results are represented in Section 5. Finally, Section 6 contains the conclusions.
2. Maximum Likelihood Estimation
Maximum-likelihood estimation (MLE) is a method of estimating the parameters of a statistical model. In general, for a fixed set of data and underlying statistical model, the method of maximum likelihood selects values of the model parameters that produce a distribution that gives the observed data the greatest probability (i.e., parameters that maximize the likelihood function). In essence the method selects a set of model parameters that predicts that events that occur often in the data are very likely to occur and events that occur seldom in the data are predicted to occur with small probability. Maximum-likelihood estimation gives a unified approach to estimation, which is well-defined in the case of the normal distribution and many other problems.
Suppose that is a random variable with probability distribution where is unknown parameters vector. Let be the observed values in a random sample of size . Then, the likelihood function of the sample is Note that the likelihood function is now a function of only the unknown parameters . The maximum likelihood estimators of are the values that maximize the likelihood function . In practice, it is often more convenient to work with the logarithm of the likelihood function: The maximum likelihood estimators would be found by equating the partial derivatives to zero and solving the resulting system of equations: While the method of maximum likelihood is an excellent technique for many models and a maximum likelihood estimator can be found as an explicit function of the observed data , sometimes complications arise in its use. For example, it is not always easy to maximize the likelihood function because the equation(s) obtained from partial derivatives may be difficult to solve or no closed-form solution to the maximization problem is known or available. Therefore, an MLE has to be found numerically using optimization methods.
3. General Iterated Local Search
Iterated local search is a simple but powerful metaheuristic algorithm [15]. It applies local search to an initial solution until it finds a local optimum; then it perturbs the solution and restarts local search. The importance of the perturbation is obvious: a too small perturbation might not enable the system to escape from the basin of attraction of the local optimum just found. On the other side, a too strong perturbation would make the algorithm similar to a random restart local search.
A local search is effective if it is able to find good local optima, that is, if it can find the basin of attraction of those states. When the search space is wide and/or when the basin of attraction of good local optima is small, a simple multistart algorithm is almost useless. An effective search could be designed as a trajectory only in the set of local optima , instead of in the set s of all the states.
The requirement on the perturbation of is to produce a starting point for local search such that a local optimum different from is reached. However, this new local optimum should be closer to than a local optimum produced by a random restart. The acceptance criterion acts as a counter balance, as it filters and gives feedback to the perturbation action, depending on the characteristics of the new local optimum. A high level description of ILS steps is presented in Algorithm 1.
| ||||||||||||||||
The design of ILS algorithms has several degrees of freedom in the choice of the initial solution, perturbation, and acceptance criteria.
The construction of initial solutions should be fast, and initial solutions should be a good starting point for local search. The fastest way of producing an initial solution is to generate it at random; however, this is the easiest way for problems that are unconstrained, whilst in other cases the construction of a feasible solution requires also constraint checking.
The perturbation is usually nondeterministic in order to avoid cycling. Its most important characteristic is the strength, roughly defined as the amount of changes made on the current solution. The strength parameter, , can be either fixed or variable. In the first case, the distance between perturbation and local search is kept constant, independently of the problem size. However, a variable strength is in general more effective, since it has been experimentally found that, in most of the problems, the bigger the problem size is, the larger should be the strength.
A second choice is the mechanism to perform perturbations. This may be a random mechanism, or the perturbation may be produced by a deterministic or semideterministic method.
The third important component is the acceptance criterion. Two extreme examples consist in (1) accepting the new local optimum only in case of improvement and (2) in always accepting the new solution. In-between, there are several possibilities. In this paper, we consider an annealing scheme for acceptance criterion. Simulated annealing is one of the most novel algorithms initially presented by Kirkpatrick et al. [16]. The algorithm is based upon that of Metropolis et al. which was originally proposed as a means of finding the equilibrium configuration of a collection of atoms at a given temperature [17]. Similar to other metaheuristic algorithms, it attempts to solve hard combinatorial optimization problems through controlled randomization. Sometimes simulated annealing was applied to parameters estimation [18].
The performance of ILS algorithm depends on the definition of the several control parameters. The choice of an appropriate is crucial for the performance of the algorithm. The value of parameter decrease during the search process; thus, at the beginning of the search, diversification is high and as it gradually goes on its search path, intensification becomes intense. Hereby, with the choice of appropriate , a dynamic balance is given between diversification and intensification.
4. Applying Hybrid ILS for Estimating Gamma/Gompertz Parameters
To estimate the three parameters of Gamma/Gompertz distribution, let be a random sample from the Gamma/Gompertz distribution. The log of the likelihood function is When the derivatives are equated to zero, we obtain the following equations that must be solved to find the maximum likelihood estimators of , , and: There is no closed form solution to these equations and it is very difficult to solve these equations using ordinary optimization techniques. To estimate the parameters , , and we are to maximize (or ), using hybrid iterated local search.
For a given values of and , (11) gives an expression of as a function of only which is called reduced . Although is a nonconcave function, it is obvious from expression (11) that reduced is a concave function because second partial derivative of reduced with respect to is nonpositive. Therefore, taking the first derivative of reduced with respect to and setting yield Therefore, problem reduces to obtain parameters and ; is computable from relation (12). Of course, only positive values for are acceptable. To challenge this very difficult problem, we propose a hybrid iterated local search approach. In this regard, the steps of this algorithm are briefly presented in Algorithm 2, where the following notation is used: : Initial solution : Current solution : Perturbed solution : Local optimum solution : Best solution : Value of the objective function at solution : Neighborhood step-length parameter : Perturbation step-length parameter : Geometric cooling factor : Initial temperature.The algorithm starts with an initial random solution for the problem (positive values for and ), and by initializing the so-called perturbation step-length parameter , the neighborhood step-length parameter , cooling factor and initial temperature . At each solution, parameter is calculated by relation (12) and only solutions are considered feasible for which is positive. Two-neighborhood solution is generated by right shift and left shift of current solution by an amount of , for or , at random. Feasible improving direction is detected. The generated feasible solution in improving direction replaces the current one. Generating feasible neighborhood solution at same direction continues until a local optimum solution is reached. This first local optimum solution sets as best solution . The cycle of algorithm repeats while the acceptance criterion is satisfied. At each repetition of the algorithm, a perturbed solution of the current local optimum solution is generated as follows: with generating two binary random numbers and , we select a direction for perturbation (right or left) for each of the parameters and . Perturbed solution is obtained by adding to or subtracting from current local optimum solution, , a dynamic amount, depending on the perturbation directions and , where is the perturbation step-length parameter playing an important role in our algorithm. The generated feasible solution replaces the current one. Local search procedure is applied to the newly chosen solution. We suppose that after the local optimum is reached, it is always acceptable (with checking positivity of ). After generating newly local optimum solution, , a kind of annealing schedule is considered as acceptance criterion: accept all the improving new local optima and accept also the nonimproving ones with a probability that is a function of the temperature and the difference of objective values, in formulas: It is obvious that control parameter, , is chosen with respect to the specific problem at hand. When adapting this general algorithm to a specific problem, the procedure to generate both initial and perturbed solutions is very important in addition to the control parameter. Hereby, the most important feature of this algorithm, as a metaheuristic, is the possibility of accepting a worst solution, which can allow it to prevent falling into local optimum trap. Obviously, the probability of accepting a worse solution decreases as the temperature decreases in each outer cycle. The initial temperature () should be high enough that, in the first iteration of the algorithm, the probability of accepting a worse solution is, at least, of 80% [16]. The most commonly used temperature reducing function is geometric; that is, in which is constant. The termination condition happens when the system has reached to a desired energy level.
| ||||||||||||||||||||||
5. Numerical Examples
In this section, we solve two problem instances for illustrative and comparative purposes. In addition, to know the effects of the sample size, different sample sizes have been chosen and used to solve estimate parameters. We consider two Gamma/Gompertz distributions with parameters , , and (Example 1) and , , and (Example 2), respectively. Let be the observed values from each distribution in a random sample of size , 200, and 300. The coding has been done using MatLab 7.12 and run on an Intel Core i5 CPU processor with 4 GB of RAM.
The performance of metaheuristic algorithms depends excessively on the value of their parameters. To know the effects of the perturbation rule and acceptance criterion, different approaches have been chosen and used to estimate parameters. We obtained best results with perturbation rule according to (in which is iteration counter) and acceptance criterion based on .
In this paper, the proposed algorithm parameter values are selected through the computational experiments. The experiments indicated that the best values for the parameters of the algorithms in our examples are as follows: and . is selected for the examples.
Tables 1 and 2 show the results. Figures 2 and 3 show the likelihood function for all the values evaluated during the process of maximization. The experimental results demonstrate that calibration provides high quality solutions. It is quite clear that as the sample size increases the estimation will be better. It is straightforward from the primary estimation theory that the bigger the sample size the better the estimation. With comparing the CPU time requirements, it is clear that as the sample size increases more run time is needed. It is because as the sample size increases the more complicated will be the likelihood function to maximize.
(a)
(b)
(c)
(a)
(b)
(c)
6. Conclusions
We considered the Gamma/Gompertz distribution, which is often encountered in customer lifetime and mortality risks distribution literature. We developed hybrid iterated local search to estimation of Gamma/Gompertz parameters. To improve the efficiency of the proposed algorithm, its parameters were selected through the computational experiments. The results of the computational experiment showed that calibration provides high quality solutions. The large size samples perform better than the small size ones, measured by the solution accuracy index but with more CPU time.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.