An Iterated Local Search Algorithm for Estimating the Parameters of the Gamma / Gompertz Distribution

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.


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  () =     exp (−  ) ; 0 ≤  < ∞,  > 0,  > 0, 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    () = (1/!)∫ ∞ 1 (ln )   −  −  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  (mean = /), 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  = 1, 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 () = (1/)[/( − 1)] ln  for ,  > 0,  ̸ = 1 and () = 1/ for ,  > 0,  = 1.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.

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 ⃗  = ( 1 ,  2 ,  3 , . . .,   ) is unknown parameters vector.Let  1 ,  2 , . . .,   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 ⃗  = ( 1 ,  2 ,  3 , . . .,   ).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  1 , . . .,   , 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.

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.

Applying Hybrid ILS for Estimating Gamma/Gompertz Parameters
To estimate the three parameters of Gamma/Gompertz distribution, let  1 , . . .,   be a random sample from the 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 Ln (, , )), using hybrid iterated local search.
For a given values of  and , (11) gives an expression of Ln (, , ) as a function of only  which is called reduced Ln (, , ).Although Ln (, , ) is a nonconcave function, it is obvious from expression (11) that reduced Ln (, , ) is a concave function because second partial derivative of reduced Ln (, , ) with respect to  is nonpositive.Therefore, taking the first derivative of reduced Ln (, , ) with respect to  and setting (/)  Ln (, , ) = 0 yield Therefore, problem reduces to obtain parameters  and ;  is computable from relation (12)    →   * , 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 ( 0 ) 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,   =  −1 in which  < 1 is constant.The termination condition happens when the system has reached to a desired energy level.
) s 0 ← Generate Initial Solution (s)  * ← Apply Local Search (s 0 ) While termination condition not met do   ← Apply perturbation ( * )   * ← Apply Local Search (  )  * ← Apply acceptance criterion ( * ,   * ) (12) 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: , cooling factor  and initial temperature  0 .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