Print["Obtain the MLE for GED"]
Print["    "]
Print["Mohammad Z. Raqab"]
Print["Department of Mathematics, University of Jordan"]
Print["July, 2013"]
Print["    "]
Print[alph,"   ",n,"      ",PC,"       ",PCK,"      ",PCF,"      ",PCAV,"       ",PCAV1]
Print["-----------------------------------------------------------------------------"]
lam=1.0;alph=1.0;It=5000;
Do[
Icount=0;Jcount=0;I1count=0;J1count=0;Kcount=0;K1count=0;Mcount=0;M1count=0;
Lcount=0;L1count=0;
alph1=1/alph;
a1 = PolyGamma[2];
a2 = PolyGamma[1, 1];
Do[
dist=UniformDistribution[{0,1}];
U=Table[Random[dist],{n}];
X=Sqrt[(-1)*Log[1-(U^alph1)]]/lam;
OX=Sort[X];
Fn=Table[N[i/n],{i,1,n}];
FG[x_,a_,b_]:=(1-Exp[-(b*x)^2])^a;
Fw[x_,c_,d_]:=1-Exp[-(d*x)^c];
Xbar=Sum[X[[i]]^2,{i,1,n}]/n;
E1[b_]:=-n/Sum[Log[1-Exp[-(b*X[[i]])^2]],{i,1,n}];
E2[b_]:=Sum[(X[[i]]^2)*Exp[-(b*X[[i]])^2]/(1-Exp[-(b*X[[i]])^2]),{i,1,n}];
f[b_]:=((E1[b]-1)*(b^2))*E2[b]-(n*(b^2)*Xbar)+n;
Lf[a_,b_]:=(n*Log[2])+(n*Log[a])+(2*n*Log[b])+Sum[Log[X[[i]]],{i,1,n}]-((b^2)*Sum[X[[i]]^2,{i,1,n}])+(a-1)*(Sum[Log[1-Exp[-(b*X[[i]])^2]],{i,1,n}]);
E3[d_]:=(n/Sum[X[[i]]^d,{i,1,n}])^(1/d);
w[d_]:=(n/d)+Sum[Log[X[[i]]],{i,1,n}]-((E3[d]^d)*Sum[(X[[i]]^d)*Log[X[[i]]],{i,1,n}]);
Lg[c_,d_]:=(n*Log[d])+(n*d*Log[c])+((d-1)*Sum[Log[X[[i]]],{i,1,n}])-((c^d)*Sum[X[[i]]^d,{i,1,n}]);
b0[j]=1.5;
b1[j]=b0[j]-(f[b0[j]]/f'[b0[j]]);
k=0;
While[Abs[b1[j]-b0[j]]>0.0005,
b0[j]=b1[j];
b1[j]=b0[j]-(f[b0[j]]/f'[b0[j]]);
k=k+1;
];
d0[j]=1.9;
d1[j]=d0[j]-(w[d0[j]]/w'[d0[j]]);
k=0;
While[Abs[d1[j]-d0[j]]>0.0005,
d0[j]=d1[j];
d1[j]=d0[j]-(w[d0[j]]/w'[d0[j]]);
k=k+1;
];
MLEalph[j]=E1[b1[j]];
MLEthet[j]=E3[d1[j]];

f11[j]= (a2 + (a1^2))/(d1[j]^2);
f22[j]= (d1[j]^2)/(MLEthet[j]^2);
f12[j]= (1 + PolyGamma[1])/MLEthet[j];
FIW[j]= {{f11[j], f12[j]}, {f12[j], f22[j]}};
FWe[j]= f11[j]+f22[j];

eta[y_,j]:= 1 + Log[1 - y] - (MLEalph[j] - 1)*(((1/y) - 1)*Log[1 - y]);
tau[y_,j]:= 1 + (MLEalph[j]*Log[y]);
g11[j]= 1/(MLEalph[j]^2);
g12[j] = (2/b1[j])*NIntegrate[tau[y,j]*eta[y,j]*(y^(MLEalph[j] - 1)), {y, 0, 1},
MaxRecursion->40];
g22[j] = ((4*MLEalph[j])/(b1[j]^2))*NIntegrate[(eta[y,j]^2)*(y^(MLEalph[j] - 1)),
{y, 0, 1},MaxRecursion->40];
FIG[j]= {{g11[j], g12[j]}, {g12[j], g22[j]}};
FGR[j] = g11[j]+g22[j];

FGR11[p_,j] := (-1/(2*b1[j]*(MLEalph[j]^2)))*Log[p]*(p^(1/MLEalph[j]))*
((-Log[1 - p^(1/MLEalph[j])])^(-0.5))/(1 - p^(1/MLEalph[j]));
FGR22[p_,j]:= (-1/(b1[j]^2))*((-Log[1 - p^(1/MLEalph[j])])^0.5);
FWE11[p_,j]:= (-1/(MLEthet[j]*(d1[j]^2)))*((-Log[1 - p])^(1/d1[j]))
*Log[-Log[1 - p]];
FWE22[p_,j]:= (-1/(MLEthet[j]^2))*((-Log[1 - p])^(1/d1[j]));
PTHG[p_,j]:= {{FGR11[p,j], FGR22[p,j]}};
PTHW[p_,j]:= {{FWE11[p,j], FWE22[p,j]}};
PTHGT[p_,j]:= Transpose[{{FGR11[p,j], FGR22[p,j]}}];
PTHWT[p_,j]:= Transpose[{{FWE11[p,j], FWE22[p,j]}}];
VarGR[p_,j]:= PTHG[p,j].Inverse[FIG[j]].PTHGT[p,j];
VarWe[p_,j]:= PTHW[p,j].Inverse[FIW[j]].PTHWT[p,j];
AVG[j] = NIntegrate[VarGR[p,j][[1,1]], {p, 0, 1},MaxRecursion->40];
AVW[j] = NIntegrate[VarWe[p,j][[1,1]], {p, 0, 1},MaxRecursion->40];

T[j]=Lf[MLEalph[j],b1[j]]-Lg[MLEthet[j],d1[j]];
FDiff[j]=FGR[j]-FWe[j];
VarDiff[j]=AVG[j]-AVW[j];
VarDiff1[j]=VarGR[0.5,j][[1,1]]-VarWe[0.5,j][[1,1]];

FtG=Table[FG[OX[[i]],MLEalph[j],b1[j]],{i,1,n}];
Ftw=Table[Fw[OX[[i]],d1[j],MLEthet[j]],{i,1,n}];
DG=Table[Abs[Fn[[i]]-FtG[[i]]],{i,1,n}];
DW=Table[Abs[Fn[[i]]-Ftw[[i]]],{i,1,n}];
KDG[j]=Max[DG];KDW[j]=Max[DW];
KDiff[j]=KDG[j]-KDW[j];

If[T[j]>0,Icount=Icount+1,Jcount=Jcount+1];
If[KDiff[j]<0,I1count=I1count+1,J1count=J1count+1];
If[FDiff[j]>0,Kcount=Kcount+1,K1count=K1count+1];
If[VarDiff[j]<0,Mcount=Mcount+1,M1count=M1count+1];
If[VarDiff1[j]<0,Lcount=Lcount+1,L1count=L1count+1],
{j,1,It}];
PC=N[Icount/It];
PCK=N[I1count/It];
PCF=N[Kcount/It];
PCAV=N[Mcount/It];
PCAV1=N[Lcount/It];
Print[alph,"    ",n,"     ",PC,"     ",PCK,"     ",PCF,"   ",PCAV,"     ",PCAV1],
{n,20,50,10}];