| rextreme=function(n,mu,sig)mu+sig(log(-log(1-runif(n)))) |
| dextreme=function(x,mu,sig)(1/sig)exp((x-mu)/sig-exp((x-mu)/sig)) |
| pextreme=function(x,mu,sig)1-exp(-exp((x-mu)/sig)) |
| n=13 |
| r=10 |
| data=log(c(0.22, 0.50, 0.88, 1.00, 1.32, 1.33, 1.54, 1.76, 2.50, 3.00)) |
| data1=sort(data)[1:r] |
| Xr=data1[r] |
| # AMLE of normal |
| pr=r/(n+1) |
| qr=1-pr |
| invpr=qnorm(pr,0,1) |
| alpha=dnorm(invpr)((1+invpr2)qr-invprdnorm(invpr))/(qr2) |
| beta=dnorm(invpr)(dnorm(invpr)-invprqr)/(qr2) |
| A=sum(data1)+beta(n-r)Xr |
| M=r+beta(n-r) |
| C=(n-r)alpha |
| D=AC/M-CXr |
| E=sum(data12)+(n-r)beta(Xr2)-A2/M |
| sigma_hat=(-D+(D2+4rE)(1/2))/(2r) |
| u_hat=A/M+Csigma_hat/M |
| ## MLE of normal |
| L2=function(x) |
| u=x[] |
| sigma=x[] |
| -(prod(dnorm(data1,u,sigma))(1-pnorm(Xr,u,sigma))(n-r)) |
| Ans1=optim(c(u_hat,sigma_hat),L2) |
| uN=Ans1$par[] |
| sigmaN=Ans1$par[] |
| ## AMLE of extreme value |
| prE=r/(n+1) |
| qrE=1-prE |
| alphar=1+log(qrE)-log(qrE)log(-log(qrE)) |
| betar=-log(qrE) |
| SUMbeta=SUMbetaX=SUMalpha=SUMalphaX=SUMbetaX2=0 |
| for(h in 1:r) |
| pi=h/(n+1) |
| qi=1-pi |
| alphai=1+log(qi)-log(qi)log(-log(qi)) |
| betai=-log(qi) |
| SUMbeta=SUMbeta+betai |
| SUMbetaX=SUMbetaX+betaidata1[h] |
| SUMalpha=SUMalpha+alphai |
| SUMalphaX=SUMalphaX+alphaidata1[h] |
| SUMbetaX2=SUMbetaX2+betai((data1[h])2) |
| M=SUMbeta+betar(n-r) |
| B=(SUMbetaX+(n-r)betarXr)/M |
| C=(SUMalpha-(n-r)(1-alphar))/M |
| D=-(n-r)Xr+(n-r)alpharXr+SUMalphaX-BCM |
| E=(n-r)betar(Xr2)+SUMbetaX2-M(B2) |
| sigma_hat=(-D+(D2+4rE)(1/2))/(2r) |
| u_hat=B-Csigma_hat |
| ## MLE of extreme value |
| L9=function(x) |
| u=x[] |
| sigma=x[] |
| -(prod(dextreme(data1,u,sigma))(1-pextreme(Xr,u,sigma))(n-r)) |
| Ans9=optim(c(u_hat,sigma_hat),L9) |
| uE=Ans9$par[] |
| sigmaE=Ans9$par[] |
| ## Model Selection Approaches |
| Dsp1=Dsp2=D1=D2=array() |
| for(j in 1:r) |
| L7=factorial(n)/(factorial(n-r)) |
| LN=L7prod(dnorm(data1,uN,sigmaN))(1-pnorm(Xr,uN,sigmaN))(n-r) |
| LE=L7prod(dextreme(data1,uE,sigmaE))(1-pextreme(Xr,uE,sigmaE))(n-r) |
| Dsp1[j]=(2/pi)abs(asin(sqrt((j-0.5)/n))-asin(sqrt(pnorm(data1[j],uN,sigmaN)))) |
| Dsp2[j]=(2/pi)abs(asin(sqrt((j-0.5)/n))-asin(sqrt(pextreme(data1[j],uE,sigmaE)))) |
| D1[j]=(2/pi)abs(((j-0.5)/n)-pnorm(data1[j],uN,sigmaN))+0.5/n |
| D2[j]=(2/pi)abs(((j-0.5)/n)-pextreme(data1[j],uE,sigmaE))+0.5/n |
| DspN=max(Dsp1) |
| DspE=max(Dsp2) |
| #DN=max(D1) |
| #DE=max(D2) |
| ## The Taylor series prediction |
| AMLP_E=function(r,n) |
| pr=r/(n+1) |
| qr=1-pr |
| Xs2=array() |
| for(i in 1:(n-r)) |
| s=r+i |
| ps=s/(n+1) |
| qs=1-ps |
| alphas=1+log(qs)-log(qs)log(-log(qs)) |
| betas=log(qs) |
| gamma1=qslog(qs)((qr-qs)(-1+(1+log(qs))log(-log(qs)))+qslog(qs)log(-log(qs))- |
| qrlog(qr)log(-log(qr)))/((qr-qs)2) |
| rou1=qslog(qs)(-(1+log(qs))(qr-qs)-qslog(qs))/((qr-qs)2) |
| v1=qslog(qs)qrlog(qr)/((qr-qs)2) |
| A=s-r-1 |
| B=Arou1+Av1+betas+(n-s)betas |
| C=Agamma1+alphas-(n-s)+(n-s)alphas |
| D=Arou1+betas+(n-s)betas |
| Xs2[i]=-Av1data1[r]/D+uEB/D-sigmaEC/D |
| Xs2[which(Xs2<=data1[r])]=data1[r] |
| Xs2 |
| AMLP_E(r,n) |
