???? 1

Size: px

Start display at page:

Download "???? 1"

소라 정
3 years ago
Views:

1 The Korean Journal of Applied Statistics (2013) 26(5), DOI: Objective Bayesian Estimation of Two-Parameter Pareto Distribution Young Sook Son a,1 a Department of Statistics, Chonnam National University (Received July 15, 2013; Revised August 23, 2013; Accepted October 15, 2013) Abstract An objective Bayesian estimation procedure of the two-parameter Pareto distribution is presented under the reference prior and the noninformative prior. Bayesian estimators are obtained by Gibbs sampling. The steps to generate parameters in the Gibbs sampler are from the shape parameter of the gamma distribution and then the scale parameter by the adaptive rejection sampling algorism. A numerical study shows that the proposed objective Bayesian estimation outperforms other estimations in simulated bias and mean squared error. Keywords: 2-parameter Pareto distribution, L-moment estimation, maximum likelihood estimation, noninformative prior, reference prior, objective Bayesian estimation, Gibbs sampling, adaptive rejection sampling. 1. 서론 Pareto (1897) 에의해처음소개되었던파레토분포 (Pareto distribution) 는인구규모, 주가변동, 보험 손실, 보험금청구, 소득, 급여, 인터넷트래픽, 최대풍속, 최대강우량, 지진혹은산불과같은자연재 해의규모, 유전지역에서오일매장량등과같이인구, 사회, 경제, 자연환경적현상등많은분야의자 료에적용되는확률모형으로중요한역할을해왔다 (Crovella 등, 1998; Harris 등, 2000; Smith, 2003; Embrechts 등, 1997; Jackson 과 Kagan, 1999; Cumming, 2001; Arnold, 1983; Arnold 와 Press, 1989). 파레토분포는지수분포와모양은비슷하나꼬리가더두터우므로극대값의발생이상대적으로더많이 발생하는자료의적용에유용하다. 2- 모수파레토 pdf(probability density function) 는다음과같이정 의된다. f(x) = αβα, x β, α, β > 0, (1.1) xα+1 여기서 α 는분포의모양을결정하는형태 (shape) 모수이고, β 는척도 (scale) 모수이다. 분포가시작하 는위치를나타낸다는의미에서 β 를위치 (location) 모수라고도한다. 식 (1.1) 처럼정의되는파레토분 포를 Pareto(α, β) 분포라나타내자. Figure 1.1 은모수 α, β 의여러값에따른 Pareto(α, β) pdf 를보 This research was supported by Basic Science Research Program through the National Research Foundation of Korea(NRF) funded by the Ministry of Education, Science and Technology(NRF ). 1 Professor, Department of Statistics, Chonnam National university, 300 Yongbong-Dong, Buk-Gu, Kwangju , Korea. ysson@jnu.ac.kr

2 714 Young Sook Son Figure 1.1. Pareto pdf for various α and β 여준다. 형태모수 α는분포의산포 (dispersion) 를나타내는척도로서 α가작을수록분포의꼬리가두터워지며산포가커진다. 반대로 α가커지면 β에서확률 1을가지는퇴화분포 (degenerate distribution) 가된다. 파레토분포는 Pickands (1975) 가제안한일반화파레토분포의특별한형태이다. 2-모수파레토분포가하한 (lower limit) 만갖는반면에일반화파레토분포는상한 (upper limit) 도허용한다. 그동안 2- 모수파레토분포의모수추정에관한많은연구가있어왔다. Singh와 Guo (1995) 는전통적인추정방법으로적률법, 확률가중법, 최대엔트로피법, 최우추정법의성능을평가하는모의실험을수행하였다. RMSE(root mean square error) 의기준하에서적률법은다른세개의추정법보다현저하게큰 RMSE 값을주었으며소표본및대표본에서최우추정법이가장우수하였다. 일반적으로최우추정법을선호하는대표적인이유는확률분포가정규조건 (regularity condition) 을만족하는경우에최우추정량이점근정규성을가지므로대표본에서모수의함수에대한구간추정을가능하게한다. 일반적으로확률분포의지지 (support) 가모수에종속하는경우에는정규조건을만족시키지못한다 (Casella와 Berger, 2002). 파레토분포가이에해당한다. 이에비해서깁스샘플링 (Gibbs sampling) 에의한베이지안추론에서는결합사후분포로부터주변사후분포를구할때적분계산에의하지않고깁스샘플러 (Gibbs sampler) 에서모수를생성하여주변사후분포를얻을수있으며, 생성된모수들의표본으로부터모수뿐만아니라모수들의함수에대한추정량의분포도쉽게얻을수있다. 따라서모수의함수들에대한구간추정도가능하다. Arnold와 Press (1989) 는확률표본및중도절단표본의경우에공액사전분포 (conjugate prior) 의가정하에서베이지안추정을제안하였다. Berger (2006) 는베이지안추론에서사전분포에대한정보가부족할때주관적사전분포의사용은좋지않은추론결과를줄수있다고지적하였다. 사전분포에대한정보가없을때디폴트사전분포 (default prior) 라고하는무정보사전분포 (noninformative prior) 하에서수행하는객관적베이지안추론 (objective Bayesian inference) 이대안이될수있다고하였다. Fu 등 (2012) 은전진적인유형 II 중도절단자료 (progressive Type II censored data) 에대하여디폴트사전분포인준거사전분포 (reference prior) 하에서객관적베이지안추론을제안하였다. 본연구에서는 Fu 등 (2012) 이사용한준거사전분포의가정하에서 2-모수파레토분포에대한객관적베이지안모수추정절차를제안한다. 제안된베이지안추정법은모의실험과자료분석을통하여 L-적률추정법및최우추정법과비교평가된다.

3 Bayesian Estimation of Pareto Distribution 고전적추정방법본절에서는고전적추정방법으로서 L-적률추정법및최우추정법에의하여 2-모수파레토분포의모수 α 및 β에대한추정량을계산한다. X = {X 1, X 2,..., X n } 를식 (1.1) 의 Pareto(α, β) 분포로부터추출한확률표본이라하고, X (1) X (2) X (n) 를만족하는 X (j) 를 j번째순서통계량이라고하자 L- 적률추정법 Pareto(α, β) 분포의 cdf(cumulative distribution function) 는다음과같다. ( ) α β F (x) = 1, x β, α, β > 0. (2.1) x Greenwood 등 (1979) 이처음제안하여수문학분야에서널리사용되어온확률가중적률 (probability weighted moment) 추정량은식 (2.2) 로부터, 그리고 Hosking (1990) 이제안한 L- 적률추정량은식 (2.3) 으로부터얻을수있다. L- 적률은순서통계량을사용하여분포의위치, 산포, 왜도, 첨도를정의할 수있으며기존의적률에비하여분포의특징을잘나타내주어매우유용하게사용되어오고있다. ( ) s α s = E [X{1 F (X)} s s ] = ( 1) r β r, s = 0, 1, 2,..., (2.2) r r=0 β r = E [X{F (X)} r ], r = 0, 1, 2,.... (2.3) 식 (2.3) 에식 (2.1) 의 cdf 를대입한적분계산은쉽지않다. 따라서 α s 를통하여 β r 를구한다. 식 (2.2) 로 부터다음을얻는다. α s = αβ α(s + 1) 1, α > 1 s + 1, α 0 = β 0, α 1 = β 0 β 1. (2.4) Pareto(α, β) 분포는모수의개수가두개이므로다음과같이처음두개의 L- 적률 λ 1 과 λ 2 가필요하다. λ 1 = β 0, λ 2 = 2β 1 β 0. (2.5) 식 (2.4) 와식 (2.5) 의관계에서모 L- 적률 λ 1, λ 2 를표본 L- 적률 l 1, l 2 로각각대체하여얻어진 α 와 β 에 대한 L- 적률추정량 (L-moment estimator; LME) 은다음과같다. ˆα LME = l 1 + l 2 2l 2, ˆβLME = (ˆα LME 1)l 1 ˆα LME, 여기서 l 1 = 1/n n X (j), l 2 = 1/{n(n 1)} n (2j n 1)X (j) 최우추정법 Pareto(α, β) 분포의우도함수 (likelihood function), ( n ) α L(α, β X) α n β nα X j (2.6) 를최대로해주는최우추정량 (maximum likelihood estimator; MLE) 은다음과같다. n ˆα MLE = n, ˆβMLE = X (1). (2.7) ln X j n ln X (1)

4 716 Young Sook Son 3. 베이지안추정 3.1. 공액사전분포하에서베이지안추정 본절에서는 Arnold 와 Press (1983) 가제안한공액사전분포가정하에서의베이지안모수추정결과를 소개한다. 식 (3.1) 과같이정의되는 Pareto(α, β) 분포에대한공액사전분포는 Lwin (1972) 에의해처 음제안되었고 Arnold 와 Press (1983) 에의해일반화되었다. α Gamma(a 0, b 0), β 1 α Pareto(αc 0, d 0), (3.1) 여기서 α 의사전 pdf 는 π(α) α a 0 1 exp( b 0 α) 이고, a 0, b 0, c 0, d 0 는양의값을가지는초모 수 (hyperparameter) 이다. 공액사전분포의식 (3.1) 하에서주변사후분포는다음과같이얻어진다.. α X Gamma(a 0 + n, b 0 + v), β 1 α, X Pareto(α(c 0 + n), min(d 0, X (1) )), (3.2) 여기서 v = n ln Xj + c0 ln d0 (n + c0) ln(min(d0, X (1))). 식 (3.2) 로부터 α 및 β 의각사후최빈값 (posterior mode) 을공액사전분포하에서의베이지안추정량 으로하였을때다음과같이얻어진다 무정보사전분포하에서베이지안추정 ˆα conjugate = a0 + n 1, ˆβconjugate = min(d 0, X (1) ). (3.3) b 0 + v 본절에서는 Fu 등 (2012) 이제안한객관적사전분포로서다음의준거사전분포 (reference prior) 의가정 하에서객관적베이지안모수추정절차를제안한다. π(α, β) 1 αβ. (3.4) 결합사후분포는우도함수식 (2.6) 과사전분포식 (3.4) 를곱하여다음과같이얻어진다. ( n ) α π(α, β X) α n 1 β nα 1 X j. (3.5) 사전분포식 (3.4) 는부적절사전분포 (improper prior) 이다. 부적절사전분포를사용하는경우에사후 분포가 pdf 로서부적절한경우에는베이지안추정이불가능하다. 그러나다음에의해식 (3.5) 에서정 의되는사후분포는 pdf로서적절함을알수있다. x(1) [ { n }] 1 π(α, β X)dβdα = n α(n 1) 1 exp α ln x j n ln x (1) dα Γ(n 1) = { } n <. n 1 n ln x j n ln x (1) 식 (3.5) 의결합사후분포로부터완전조건부사후분포는다음과같이정의된다. [ { n }] π(α β, X) α n 1 exp α ln x j n ln β. (3.6) π(β α, X) β nα 1. (3.7)

5 Bayesian Estimation of Pareto Distribution 717 Initial step: Adopt β (0) as the initial value of β. Replication step: for i = 0, 1, 2,..., I, ( n - Generate α (i+1) Gamma n, ln x j n ln β ). (i) - Generate β (i+1) by ARS algorithm given α (i+1). Figure 3.1. Gibbs Sampler for generating α and β 기각표집법 (rejection sampling) 은 pdf 로부터확률변수를생성하는일반적인방법이다. 기각표집법을 사용하여정규화상수위까지만명시되어있는, 즉, pdf 의 kernel 만주어진 kernel pdf 로부터변수를생 성하고자하는경우에필요한 kernel pdf 의덮개함수 (envelope function) 를찾는것은현실적으로어 렵다. Gilks 와 Wild (1992) 가제안한적응기각표집법 (adaptive rejection sampling; ARS) 은 kernel pdf 의로그 - 오목성 (log-concavity) 을이용하여 kernel pdf 의로그함수에접선을연결하여만든조각직 선덮개함수 (piecewise linear envelope function) 를사용한다. 식 (3.6) 은 Gamma pdf 의커널이며, 식 (3.7) 에서 n > 1/α 인경우에 2 ln π(β α, X)/ β 2 = (nα 1)/β 2 < 0 가되어 π(β α, X) 는로그 - 오목성을갖는다. n > 1/α 는현실적으로만족되는 조건이다. 따라서깁스샘플러를 Figure 3.1 과같이구성한다. 4. 수치분석본연구에서제안한객관적베이지안모수추정법의평가를위하여모의실험과자료분석을수행하였다. 모든계산은 R package(version ) 를사용하여이루어졌다. 파레토, 감마확률변수생성을위하여 rpareto 및 rgamma 함수를각각사용하였고 ARS 알고리즘을구현하기위하여 ars.new 함수를사용하였다. 또한 HPD(highest posterior density) interval을구하기위하여 HPDinterval 함수를사용하였다. HPDinterval 함수는경험적 cdf로부터지정된확률을가지는모든가능한구간중에서가장길이가짧은구간을 HPD interval로계산한다 모의실험척도모수 β는척도불변성 (scale invariance) 을가지므로 β = 1로놓는다. 그리고 α = 0.5, 1, 3, 5, 10인 Pareto(α, β) 로부터표본의크기가 n = 15, 30, 100인자료세트를각각 1,000회반복생성하여모의실험을하였다. n > 2이므로 ARS 알고리즘의사용이가능하다. 깁스샘플링은 11,000회반복수행하여처음 1,000회를버리고나머지 10,000회생성된 α 및 β값들로사후분포를구성하였다. 시험적인모의실험에서 α 및 β의사후분포는각각대칭및좌비대칭분포의형태를가지고있었다. 따라서 α 및 β 의각사후분포의평균과최빈값을베이지안추정량으로사용하기로하였다. Table 4.1에는 L- 적률추정량 (LME), 최우추정량 (MLE), 준거사전분포가정하의객관적베이지안추정량 (BAYES) 들의 bias와 MSE(mean squared error) 가제시되어있다. 모의실험결과는다음과같이요약된다. 첫째, 모형및표본의크기에관계없이 bias와 MSE의기준하에서 BAYES, MLE, LME 순으로우수하였다. BAYES는소표본 (n = 15) 에서 MLE보다더우수하고대표본 (n = 100) 에서는비슷하다. 둘째, α의값이커질수록분포의산포가커지므로 α 추정치들의 bias와 MSE도점점커진다. 셋째, 표본의크기가커짐에따라 bias와 MSE는점점작아진다.

6 718 Young Sook Son Table 4.1. Simulated bias and MSE for five Pareto(α, β) models estimation n = 15 n = 30 n = 100 method bias MSE bias MSE bias MSE LME α = 0.5 MLE BAYES LME β = 1.0 MLE BAYES LME α = 1.0 MLE BAYES LME β = 1.0 MLE BAYES LME α = 3.0 MLE BAYES LME β = 1.0 MLE BAYES LME α = 5.0 MLE BAYES LME β = 1.0 MLE BAYES LME α = 10.0 MLE BAYES LME β = 1.0 MLE BAYES Table 4.2. Annual wage data(unit: in multiples of 100 U.S. dollars) 자료분석 Table 4.2의자료는 Arnold와 Press (1989) 가공액사전분포의가정하에서제안한베이지안추정절차에의해 Pareto(α, β) 분포에적합하였던연간임금 (annual wage) 자료이다. 이자료의확률히스토그램인 Figure 4.1을보면꼬리부분에과도하게자료가분포하고있음을볼수있다. 그들이가정한공액사전분포의초모수는식 (3.1) 에서 a 0 = 6.02, b 0 = 3.17, c 0 = , d 0 = 이다. 이렇게가정된공액사전분포하에서베이지안추정치는식 (3.3) 에의해 ˆα conjugate = , ˆβ conjugate = 101로얻어졌다.

7 Bayesian Estimation of Pareto Distribution 719 Figure 4.1. Pareto density fitted to annual wage data Figure 4.2. Posterior distribution of α and β Table 4.3. Descriptive summary for the posterior distribution of α and β parameter mean SD median mode 95% HPD interval [lower, upper] α β 본연구에서제안한객관적베이지안모수추정을위하여 Figure 3.1의 Gibbs sampler를 I = 11,000회반복수행하여처음 1,000회를버리고나머지 10,000회반복생성한 α와 β값들로사후분포를구성하였다. 초기값 β (0) 는 ˆβ LME 를사용하였다. Figure 4.2에서보여주는 α 및 β의사후분포로부터각모수의베이지안추정치로서 α에대해서는사후평균값을 β에대해서는사후최빈값을사용하기로한다. Table 4.3은 α 및 β에대한사후분포의기술통계분석결과를보여준다. 깁스샘플링의매반복에서표집된 α 및 β의시계열도표인 Figure 4.3에서이시계열들은각각세로축에수평참조선으로표시된 α의사후평균값과 β의사후최빈값으로수렴하고있음을알수있다. 또한 α 및 β의각시계열에대한제1차자기상관계수는각각 및 0.043으로깁스샘플들이상호관련성이없음을알수있다. Table 4.4는모수 α와 β에대한각추정치와각추정된모형에대한적합도검정으로서 χ 2 -test와 KStest(Kolmogrov Smirnov test) 의 p-value를제시하고있고, Figure 4.4는 Q-Q plot(quantile-quantile

720 Young Sook Son Figure 4.3. Gibbs sampling iterations plot for α and β Table 4.4. Results fitted to annual wage data estimation parameter estimate goodness of fit test(p value) method α β χ 2 -test K-S test LME 6.

8 720 Young Sook Son Figure 4.3. Gibbs sampling iterations plot for α and β Table 4.4. Results fitted to annual wage data estimation parameter estimate goodness of fit test(p value) method α β χ 2 -test K-S test LME MLE BAYES-conjugate BAYES-default Figure 4.4. Q-Q plot plot) 을그리고 Figure 4.5는 P-P plot(probability-probability plot) 을보여준다. 이상과같은모형적합결과로부터객관적베이지안추정이최우추정보다다소더적합한그러나대체로비슷한적합도를보여주고있으며, 공액사전분포하의베이지안추정은이들보다적합도가떨어짐을알수있다. 원자료의확률히스토그램과최우추정및베이지안추정에의해추정된파레토 pdf를동시에보여주고있는 Figure 4.1에서꼬리부분에과도하게분포되어있는자료는 Figure 4.4의 QQ-plot에서극대분위수들이직선에서상당히이탈한결과와일치하고있음을알수있다.

9 Bayesian Estimation of Pareto Distribution 721 Figure 4.5. P-P plot 5. 결과본연구에서는 2-모수파레토분포에대해무정보사전분포인준거사전분포의가정하에서객관적베이지안모수추정절차를제안하였다. 베이지안추정은깁스샘플링에의해서수행된다. 깁스샘플러에서형태모수는감마분포로부터생성되고, 척도모수는 Gilks와 Wild (1992) 의 ARS 알고리즘에의해생성된다. 제안된베이지안모수추정절차는모의실험과자료분석에서기존의추정방법들인 L-적률추정법, 최우추정법, 그리고공액사전분포하의주관적베이지안모수추정법과비교된다. 제안된객관적베이지안추정법이최우추정법보다다소더좋은결과를보여주었다. 본논문에모의실험결과를제시하지는않았으나주관적베이지안추정의경우에공액사전분포의초모수선택의주관적판단에따라모형적합결과는매우좌우되며오히려오차가크게발생하게하는결과를준다. References Arnold, B. (1983). Pareto Distributions, Fairland, MD: International Co-operative Publishing House. Arnold, B. C. and Press, S. J. (1983). Bayesian inference for Pareto populations, Journal of Econometrics, 21, Arnold, B. C. and Press, S. J. (1989). Bayesian estimation and prediction for Pareto data, Journal of the American Statistical Association, 84, Berger, J. O. (2006). The case for objective Bayesian analysis (with discussion), Bayesian Analysis, 1, Casella, G. and Berger, R. L. (2002). Statistical Inference, Duxbury Press. Crovella, M. E., Taqqu, M. S. and Bestavros, A. (1998). Heavy-Tailed Probability Distributions in the World Wide Web, A Practical Guide to Heavy Tails, Adler, R. J., Feldman, R. E., and Taqqu, M. S., Editors, Birkhauser, Boston, MA, Cumming, S. G. (2001). A parametric model of the fire-size distribution, Canadian Journal of Forest Research, 31, Embrechts, P., Kluppelberg, C. and Mikosch, T. (1997). Modeling Extremal Events, Springer, New York. Fu, J., Xu, A. and Tang, Y. (2012). Objective Bayesian analysis of Pareto distribution under progressive Type-II censoring, Statistics and Probability Letters, 82, Gilks, W. R. and Wild, P. (1992). Adaptive rejection sampling for Gibbs sampling, Applied Statistics, 41,

10 722 Young Sook Son Greenwood, J. A., Landwehr, J. M., Matalas, N. C. and Wallis, J. R. (1979). Probability weighted moments: Definition and relation to parameters of several distributions expressible in inverse form, Water Resources Research, 15, Harris, C. M., Brill, P. H. and Fischer, M. J. (2000). Internet-type queues with power-tailed interarrival times and computational methods for their analysis, INFORMS Journal on Computing, 12, Hosking, J. R. M. (1990). L-Moments: Analysis and estimation of distributions using linear combinations of order statistics, Journal of the Royal Statistical Society: Series B, 52, Jackson, D. D. and Kagan, Y. Y. (1999). Testable earthquake forecasts for 1999, Seismological Research Letters, 70, Johnson, N. L., Kotz, S. and Balakrishnan, N. (1995). Continuous Univariate Distributions, 2nd ed., Wiley, NY. Lwin, T. (1972). Estimation of the tail of the Paretian law, Scandinavian Actuarial Journal, 55, Pareto, V. (1897). Cours d éeconomic Politique, F.Pichou, Paris. Pickands, III, J. (1975). Statistical inference using extreme order statistics, The Annals of Statistics, 3, Singh, V. P. and Guo, H. (1995). Parameter estimations for 2-parameter Pareto distribution by Pome, Water Resources Management, 9, Smith, R. L. (2003). Extreme Values in Finance, Telecom-Munications and the Environment, Finkenstadt, B. and Rootzen, H., Editors, Chapman and Hall/CRC Press, London.

11 Bayesian Estimation of Pareto Distribution 모수파레토분포의객관적베이지안추정 손영숙 a,1 a 전남대학교통계학과 (2013 년 7 월 15 일접수, 2013 년 8 월 23 일수정, 2013 년 10 월 15 일채택 ) 요약본연구에서는 2- 모수파레토분포에대해무정보사전분포인준거사전분포의가정하에서객관적베이지안모수추정절차를제안하였다. 베이지안추정은깁스샘플링에의해서수행된다. 깁스샘플러에서모수생성하는방법은형태모수는감마분포로부터생성하고척도모수는적응기각표집알고리즘에의해생성한다. 제안된베이지안모수추정절차는모의실험과자료분석에서기존의추정방법들인 L- 적률추정법, 최우추정법, 공액사전분포하의주관적베이지안모수추정법과비교된다. 주요용어 : 2- 모수파레토분포, L- 적률추정법, 최우추정법, 무정보사전분포, 준거사전분포, 객관적베이 지안추정, 깁스샘플링, 적응기각표집. 본연구는 2011 년도정부 ( 교육과학기술부 ) 의재원으로한국연구재단의기초연구사업지원을받아수행된것임 (NRF ). 1 ( ) 광주광역시북구용봉동 300 번지, 전남대학교, 통계학과. ysson@jnu.ac.kr

???? 1

The Korean Journal of Applied Statistics (2014) 27(1), 13 20 DOI: http://dx.doi.org/10.5351/kjas.2014.27.1.013 Maximum Tolerated Dose Estimation by Stopping Rule and SM3 Design in a Phase I Clinical Trial