Parameter estimation of Gumbel distribution using Quasi-Newton Broyden Fletcher Goldfarb Shanno (BFGS) method and its application on data of daily precipitation in Purworejo regency

Extreme data on an observation can occur due to rare events in the observation, and therefore they should be examined. One of methods to collect extreme data is block maxima. The distribution of extreme datasets collected using block maxima method is called extreme value distribution. Gumbel distribution is defined as extreme value distribution with two parameters. It is difficult to determine exact values in the parameter estimation of Gumbel distribution using maximum likelihood (ML) method. Therefore, the present research seeks to estimate the parameters using approximate solution of BFGS quasi-Newton method. The BFGS method is one of the most popular members of quasi Newton method. The purpose of this research was to determine the parameter estimation of Gumbel distribution with quasi Newton BFGS method. The quasi Newton BFGS method is a numerical method used for nonlinear function optimization without constraint so that the method can be used for parameter estimation from Gumbel distribution whose distribution function is in the form of exponential double function. The parameter estimation of the Gumbel distribution by numerical approach using the quasi Newton BFGS method is done by calculating the parameter values that make the distribution function maximum. This research is a theory research and application by studying several journals and textbooks. The results of this research we obtained the quasi Newton BFGS algorithm and estimation of Gumbel distribution parameters. Such method was applied to data of daily precipitation in Purworejo regency for the purpose of estimating the distribution parameters. The parameter estimation using ML and the quasi-Newton BFGS results in β^=βk+dkgk and γ^=γk+dkgk the Gumbel distribution for period I, β^=63.436 and γ^=49.8866 and for period II, β^=53.3121 and γ^=48.9705. This indicates that high intensity and range of precipitation have decreased.


Introduction
In observation, several data may sometimes appear to be different from or inconsistent with the rest of the data. These data are termed outliers ( [1]). They are defined as an observation (or subset of observations) which does not (or do not) follow the pattern of the main bulk of data, or which deviates (or deviate) so much from primary data ( [2]). If outliers have a greater value than other data, they belong to extreme data. The extreme data will influence the analysis and conclusion drawing. Extreme data occur due to either observation or measurement error, and unpredicted events in analysis and  [3][4]. They can provide important information which other data cannot give, and therefore their existence should be further examined ( [3]).
One method to collect extreme data by taking the maximum value in a given period is block maxima ( [4][5][6]). The given period is called block. It can be weekly, monthly, bimonthly, trimonthly, and so on. The distribution of extreme datasets collected using block maxima method is termed extreme value distribution. The Gumbel distribution is the extreme value distribution with two parameters ( [7]). The parameter estimation of the Gumbel distribution is required to analyze extreme data.
A research on the parameter estimation of the Gumbel distribution was conducted by [8,9]. Reference [9] compared method of moments, probability weighted moments (PWM), principle of maximum entropy (POME), and ML to decide the best method for the parameter estimation of the Gumbel distribution with such criteria as bias, mean-square error (MSE), and efficiency. Such research indicated that ML is the best method due to its efficiency and its smallest MSE.
Reference [8] compared three methods to estimate parameters of Gumbel distribution: method of moments, ML, and PWM. It is in fact difficult to determine exact values using ML and method of moments, and therefore numerical approximation is required to estimate the parameters. The research, however, did not mention the numerical approximation used. The research demonstrated that both PWM and ML are accurate for all sample sizes. The aforementioned researches proved that ML is an appropriate method to estimate parameters of Gumbel distribution. The parameter estimation using ML, however, presents a drawback that finding exact solution is difficult. For that reason, numerical approximation is highly required.
The BFGS quasi-Newton method is presented for approximate solution, which is so-called numerical method. It is used for unconstrained optimization of nonlinear function, and therefore it can be applied to parameters of Gumbel distribution. The BFGS quasi-Newton method is developed from Newton's method. According to [10], Newton's method makes use of the second derivative to compute the change of parameter value in each iteration. This method is later modified with the addition of step length to ensure the convergence in case that the second derivative requires complex computation. In short, Newton's method is modified by updating the second derivative in each iteration. The objective of this research is to estimate parameters of Gumbel distribution using ML and the BFGS quasi-Newton methods.

Research Method
The present research was conducted by examining the BFGS quasi-Newton method used to estimate parameters of Gumbel distribution. Such method was later applied to data of precipitation in Purworejo regency. Data were obtained from Purworejo Department of Energy and Mineral Resources (Badan Energi dan Sumber Daya Mineral-ESDM). Below are procedures carried out in the research: (1) Estimate the Gumbel distribution parameters by the MLE method. Following steps taken to estimate distribution parameters Gumbel with MLE method. a. Determining the likelihood function of the Gumbel distribution pdf as in equation (2). b. Calculates ln of the likelihood function. c. Determining the first derivative of ln -likelihood against each parameters that are estimated and maximize. (2) In step (1) we get parameter estimation and that it is a system of nonlinear equations which is difficult to determine its exact value a numerical approximation with quasi Newton BFGS is required.

Parameter Estimation
Let = { 1 , 2 , … , } denotes extreme datasets resulted from block maxima. The data are assumed to be Gumbel-distributed with such parameters as and . The parameters can be estimated using ML method. The parameter estimation using ML method is obtained by maximizing likelihood function. The likelihood function involves joint probability density function. The Gumbel probability density function is employed in the Gumbel distribution ([11]). It is expressed as: According to Bain and Engelhardt [14], a likelihood function is, therefore, obtained: Since the likelihood function belongs to exponential function, to simplify the computation it can be transformed into ln likelihood function, as expressed below: The ln likelihood function in equation (1) is included as maximum ln likelihood function if it fulfills ( , ) = 0 and ∂ln ( , ) = 0: Estimation of parameters and is the solution of system of equation (2). The system of equation (2) is a system of nonlinear equation of which exact solution is difficult to find. The parameter estimation of the Gumbel distribution, therefore, is applied using approximate solution of the BFGS quasi-Newton.
The parameter estimation of Gumbel distribution with BFGS quasi-Newton is carried out by determining parameter values that contribute to maximum distribution function. The gradient vector of ln likelihood function and Hessian matrix are required ( [12]). The gradient vector of ln likelihood function, ( , ) = c k , is expressed: The first element of is the first derivative of ln likelihood function with respect to parameter , while the second element is the first derivative of ln likelihood function with respect to parameter .
Hessian matrix, , is defined as a matrix of which elements involve the second derivative of ln likelihood function, and is denoted: , where The first element, 11 , is the second derivative with respect to parameter . The second element, 12 , is the partial derivative with respects to parameters and . The third element, 21 , is the partial derivative with respects to parameters and . The forth element, 22 , is the second derivative with respect to parameter . According to [13], the BFGS quasi-Newton method requires a negative definite matrix H to maximize a function by involving principal minor determinant. The Hessian matrix is negative definite if the principal minor determinant is negative. The Hessian matrix on matrix (4) is 2x2 sized matrix, and therefore the principal minors are 1 = 11 and 2 = . The determinant of 1 and 2 is then defined: Value of parameter > 0 and of parameter > 0, value of | 1 | < 0 , but value of | 2 | is difficult to mathematically determine. For that reason, Hessian matrix is negative definite if | 2 | < 0 so that ln likelihood function is defined as maximum function.

BFGS Quasi-Newton Algorithm
The following is BFGS quasi-Newton algorithm, an algorithm derived from [14], starting from the k th iteration: (1) Determining the initial value of gradient vector 0 , hessian matrix 0 , and matrix 0 . Gradient vector 0 is column vector (3), where k = 0 and Hessian matrix 0 is column vector (4), while matrix 0 is a matrix of which elements are the initial values of parameters and . stop iteration with the criterion of ‖ ‖ < error tolerance. In this case, the error tolerance provided is 0.0001. When the criterion is met, the iteration stops and the scores of parameter estimation are ̂= +1 and ̂= +1 . On the contrary, if ‖ ‖ > error tolerance, then to obtain parameter estimation, iteration is repeated from step (2) to step (8). Hence, parameter estimation is obtained, ̂= +1 and ̂= +1 in which +1 is calculated based on equation (6) and +1 is calculated based on equation (7).

Implementation
This part shows the implementation of parameter estimation on daily precipitation in Purworejo regency for Purwodadi precipitation station in the periods of January 1994 to December 2013. Extreme data were taken using block maxima with three monthly blocks. Data on daily precipitation in Purworejo Regency were divided into 2 periods of 10 years, namely Period I (January 1994 -December 2003) and Period II (January 2004-December 2013). Each period was divided into 4 trimonthly blocks. The blocks comprised January-February-March (JFM), April-May-June (AMJ), July-August-September (JAS), and October-November-December (OND). Afterwards, the maximum precipitation of each block was determined. The maximum precipitation of each block served as extreme datum and therefore for period of 1994 -2013, 40 data were obtained for Period I and 40 data were also obtained for Period II. Figure 1 shows an illustration of daily rainfall sampling with the maxima block method in Purworejo District in January period until December 2013. To see the pattern of extreme data distribution in period I and II is used histogram. Figure 2 shows that the pattern of extreme data distribution for periods I and II are different. However, the pattern of distribution can not be determined characteristics of extreme value of rainfall in Purworejo District. For it is estimated the distribution parameters to be determined values parameters.
(a) (b) Figure 2. Histogram of rainfall (a) period I and (b) period II A goodness of fit test for distribution was used to determine whether extreme data taken were Gumbel-distributed. This test was carried out in each period using Anderson-Darling test. The following was the goodness of fit test for period I and period II with the confidence level of 95%.
(1) Hypothesis, 0 : extreme data in period I and period II were Gumbel-distributed and 1 : extreme data in period I and period II were not Gumbel-distributed. (2) Rejection region: 0 was rejected if p−value < α.
(3) Test statistics: For period I, p−value = 0.250 > = 0.05 and for period II p−value = 0.120 > = 0.05. (4) Conclusion. Since p−value for both period was greater than the value of , 0 was rejected, meaning that that extreme data in period I and II were Gumbel-distributed. On the basis of the results of hypothesis test of extreme data of precipitation in Purworejo regency in period I and II which had followed Gumbel distribution, parameter estimation of Gumbel distribution was made using MLE method with approximate solution. Based on the BFGS quasi-Newton algorithm, initial value of H matrix was negative identity matrix, −I, and the error tolerance was 0.0001 for the computation of both periods. In period I, the initial values for parameter 0 = 63 and  Table 1. There are two parameters that are estimated, namely location and scale parameters. The location parameter indicates the location of the data center and the parameters scale shows the amount of data diversity ( [2]). In period I, the value its location parameter is 63.436 and its scale is 49.8866. This means in the first period of rainfall that occurred in Purworejo District on average amounted to 63.436 and the diversity of rainfall that occurred amounted to 49.8866. This indicating a lot of high rainfall is happening. In the second period the location parameter value is 53.3121 and the scale is 48.9705. This means in the second period of rainfall that occurred in the District Purworejo an average of 53.3121 and the diversity of rainfall that occurred amounted to 48.9705. In the second interval period also there are still many indications happen high rainfall. Average and rainfall diversity occurring on the second period decreased in value compared to average and rainfall diversity in period I or the pattern changes during the interval period I and II. According to the IPCC in 2001 [14] the pattern of rainfall change is reviewed from the value of concentration and diversity of rainfall data is determined by three criteria following. a. Changes in the value of centralization and data diversity. Improved concentration values data followed by an increase in the value of the resulting data diversity increasing the extreme value limits of rainfall with value limits extreme bottom of fixed rainfall. In addition, the probability of occurrence of bulk Increased rain increases and vice versa for decreasing concentration values and multiple data simultaneously. b. Changes in the value of data centering without being followed by changes in the value of diversity data. Improved data center values not followed by upgrades the value of data diversity shows an increased probability of occurrence high rainfall and decreased probability of bulk occurrence low rain without any change in the extremes of lower extremities nor above and vice versa. c. Changes in the value of data diversity without the change of concentration values followed data.
Increase in the value of data diversity not followed by upgrades the value of data concentration indicates a decrease in probability of occurrence high rainfall without any change in extreme value limits down or up and vice versa. Based on these criteria, the pattern of rainfall change in Purworejo corresponds to the first criterion. First criterion change shown in Figure 3. The parameter value shows the total average rainfall and shows the range of rainfall occurring. Decrease in value and simultaneously shows that in the Purworejo the occurrence of extreme value limits on rainfall and decline the probability of high rainfall occurrence.

Conclusion
In this paper, we estimated parameters of Gumbel distribution using ML and the BFGS quasi-Newton method were ̂= +1 and ̂= +1 . Then, Its application to real data of daily precipitation in Purworejo regency, we obtained two parameter estimation for period I, ̂= 63.436 and ̂= 49.8866, and for period II, ̂= 53.3121 and ̂= 48.9705. Location parameter, , shows the location of data center, while scale parameter, , indicates the distribution of data from data center. The decreasing values of location and scale parameters from period I to period II signify that data center and distribution in both periods have declined. The decreasing data distribution indicate that precipitation in period II is more homogenous than that in period I.