Thank you Rick your comments were so great and I applied them right away. I picked the values of (alpha1, alpha2 and C) that provide me with the highest log-likelihood. But my main concern is still there. That is getting negative values inside the square root or some times (depends on what initial I am using) the argument of the Log function is negative. causing my program to act weird and the estimates to go wild. I thought at the beginning that maybe I miss picked the initials so I wrote to you asking about that. But now I am not sure where is exactly the problem. My program is designed to calculate MLE, Moment and Bayes estimates using different loss functions. For Bayes method I am using Lindley's as approximation method. Thank you, Salah * Simulation from Bivariate Lomax;
********************************************************************************;
*Remarks: (1) X and Y are associated as long as alpha1*alpha2>0
(2) c>2;
* In order to figure out what initial values I need for the priors, I will try to use Shoukri's article page 35 and work it backword;
* to come up with a1 and a2 I will use gamma=a1/(a1+a2)so given a1 and gamma I will find a2 then I will use a1 and a2 and lambda for the simulation;
* a1 will be simulated from Gamma(a1, lambda) a2~Gamma(a2, lambda)....S=alpha1+alpha2 hence S~Gamma(a1+a2, lambda), gamma~Beta(a1,a2);
* Last updates Feb 9, 2018 ;
** Update 11 Feb 2018: ;
** In case of getting this message;
** ERROR: NEWRAP Optimization cannot be completed.;
** WARNING: Optimization routine cannot improve the function value.;
** Rick_SAS gave the following answer "Personally, I wouldn't worry about it. It is probably the result of degenerate (simulated) data that doesn't look very much like ;
** a bivariate Lomax (alpha1,alpha2,c) distribution. If you are curious, you can try increasing the number of Newton iterations and see if the problem goes away.;
** For reference visit this link: https://communities.sas.com/t5/SAS-IML-Software-and-Matrix/How-do-I-use-modules-within-Do-loop-groups/td-p/264123/page/2;
*Notice: When using a2=45 the linex won't work if m is small (m=10) if we increase m to 100 then we won't have this problem.
********************************************************************************;
Proc iml;
seed1=12345;* some times seed=2,0;
seed2=67801;
seed=0;
m=10;K=1000;
theta=0.1;
*make sure that mu1>1 and mu2>1;
mu1=2;mu2=4;
a1=mu1; a2=mu2;
z= 1;
mu=3;lambda=2 ; *Initial values for the gamma, I will assume a1=mu1 and a2=mu2;
C=1/theta;
********************** Initial values for r and S for the Bayes case *************************;
Prior_a1=J(10000,1,0); Prior_a2=J(10000,1,0);prior_g=J(10000,1,0);
Do I=1 to 10000;
Prior_a1[i]=rangam(seed,a1);
Prior_a2[i]=rangam(seed,a2);
prior_g[i]=Prior_a1[i]/(Prior_a1[i]+Prior_a2[i]);
end; *DO LOOP OF THE GAMMA;
alpha1=Prior_a1[:]; *true value for alpha1;
alpha2=Prior_a2[:]; *true value for alpha2;
gamma=alpha1/(alpha1+alpha2);
*gamma= prior_g[:];
***********************************************************************;
*************************************************************;
*** LOG likelihood
*************************************************************;
start LogLik(param) global (m,X,Y);
f=0;
alpha1=param[1];
alpha2=param[2];
c=param[3];
sum=log(1+alpha1*X+alpha2*Y) ;
f=f+m*log( alpha1*alpha2*c*(1+c) )-(c+2)*sum[+];
return(f);
finish;
**************************************************************;
C_MLE=J(k,1,0); Gamma_moment=J(k,1,0);gamma_mle=J(k,1,0);S_mle=J(k,1,0);Corr=J(k,1,0);
gamma_Sqr=J(k,1,0); gamma_Lin=J(k,1,0); gamma_pre=J(k,1,0);
X=J(m,1,0);Y=J(m,1,0);
Do N1=1 to K;
** Data simulation using this procedure ;
*Start Data(m,c,alpha1,alphau2,X,Y) global(seed1, seed2);
U1=J(m,1,0);U2=J(m,1,0);
do i=1 to m;
U1[i]=Uniform(seed1);
U2[i]=Uniform(seed2);
X[i]=( (1-U1[i])**(-1/c) -1 )/alpha1;*simulation is written differently;
Y[i]=(1+alpha1*X[i])/alpha2*( (1-U2[i])**(-1/(1+c)) -1 );
end;
**************************************************************************************;
/******* Compute the correlation *******;
var_X=X`*X-m*X[:]**2;
var_y=y`*y-m*y[:]**2;
Corr[N1]=(X`*Y-m*X[:]*Y[:])/sqrt( var_X*var_y);
/*************** Optimization Step to find MLEs *********************;
************* Constrain***********************/;
con={1.e-6 1.e-6 1.e-6, . . .};
opt={1}; /*1: find max*/
*x0={0,0,0};
x0=alpha1||alpha2||c; /*initial guess done by SAS*/
tc={10000 15000};
*************** NLPQN for the Lomax function *********************;
call NLPNRA(rc, param_ret, "LogLik", x0, opt, con,tc);
alpha1_mle=param_ret[1];
alpha2_mle=param_ret[2];
c_mle[N1]=param_ret[3];
*c _mle=corr;
Gamma_mle[N1]=alpha1_mle/(alpha1_mle+alpha2_mle);
S_mle[N1]=(alpha1_mle+alpha2_mle);
Gamma_moment[N1]=Y[:]/(X[:]+Y[:]);
**********************************************************************;
************************************************************************************;
************* Derivations For Lindleys method **********************;
************************************************************************************;
**L_ggg: derivative w.r.t gamma-gamma-gamma ;
**L_ggs: derivative w.r.t. gamma-gamma-S;
sum_ggg=( S_mle[N1]#(X-Y)/( 1+gamma_MLE[N1]#S_MLE[N1]#(X-Y)+S_MLE[N1]#Y ) )##3;
L_ggg=-2*m/(1-gamma_mle[N1])**3 + 2*m/gamma_mle[N1]**3 -2*(2+c_mle[N1])* sum_ggg[+];
*********************************************************************************************************;
sum_ggs1=( 2# S_mle[N1]**2#(X-Y)##2#(gamma_mle[N1]#(X-Y)+Y ))/( 1+gamma_MLE[N1]#S_MLE[N1]#(X-Y)+S_MLE[N1]#Y ) ##3 ;
sum_ggs2=( 2# S_mle[N1]#(X-Y)##2/( 1+gamma_MLE[N1]#S_MLE[N1]#(X-Y)+S_MLE[N1]#Y ) )##2;
L_ggs= -(2+c_mle[N1]) *(sum_ggs1[+]-sum_ggs2[+]) ;
**********************************************************************************************************;
* L_gss was missing here ;
sum_gss1= 2# S_mle[N1]#(X-Y)#( gamma_mle[N1]#(X-Y)+Y )##2 /( 1+gamma_MLE[N1]#S_MLE[N1]#(X-Y)+S_MLE[N1]#Y ) ##3 ;
sum_gss2=2# (X-Y) #( gamma_MLE[N1]#(X-Y)+Y )/( 1+ gamma_MLE[N1]#S_MLE[N1]#(X-Y)+S_MLE[N1]#Y )##2;
L_gss= -(2+c)*( sum_gss1[+] - sum_gss2[+] );
**************************************************************************************************;
**************************************************************************************************;
sum_sss=( (gamma_mle[N1]#(X-Y)+Y )/( 1+gamma_MLE[N1]#S_MLE[N1]#(X-Y)+S_MLE[N1]#Y ) )##3;
L_sss = 2*(1+2*m)/S_mle[N1]**3 - (2+c_mle[N1])*2*sum_sss[+];
L_gsg=L_ggs; L_sgg=L_ggs;
L_ssg= L_gss ;
**************************************************************************************************;
*** Joint Prior ***;
*** P[1] is the partial derivative of p w.r.t. gamma;
*** P[2] is the partial derivative of p w.r.t. s;
P=J(2,1,0);
P[1]=(a1-1)/gamma_mle[N1]-(a2-1)/(1-gamma_mle[N1]);
p[2]=(mu-1)/S_mle[N1]-lambda;
******** Fisher Matrix *****;
sigma= J(2,2,0);
Fisher=J(2,2,0);
Sum_I11=S_mle[N1]**2#(X-Y)##2/( 1+gamma_MLE[N1]#S_MLE[N1]#(X-Y)+S_MLE[N1]#Y )##2;
*print sigma Sum_I11;
Fisher[1,1]= -m/(1-gamma_mle[N1])**2 - m/gamma_mle[N1]**2 + (2+C_mle[N1])*Sum_I11[+];
Sum_I12= - S_mle[N1]#(X-Y)#( gamma_mle[N1]#(X-Y)+Y )/ ( 1+gamma_MLE[N1]#S_MLE[N1]#(X-Y)+S_MLE[N1]#Y )##2
+ ( X-Y )/ ( 1+gamma_MLE[N1]#S_MLE[N1]#(X-Y)+S_MLE[N1]#Y );
Fisher[1,2]= -( 2+c_mle[N1] )* Sum_I12[+];
Fisher[2,1]=Fisher[1,2];
Sum_I22= ( gamma_mle[N1]#(X-Y)+Y )##2/( 1+gamma_MLE[N1]#S_MLE[N1]#(X-Y)+S_MLE[N1]#Y )##2 ;
Fisher[2,2]= - (1+2*m)/S_MLE[N1]**2 +( 2+c_mle[N1] )* Sum_I22[+];
sigma=inv(fisher);
******************************************;
*******************************************************************************;
***************************** Square Error Loss ******************************;
*******************************************************************************;
gamma_Sqr[N1]= gamma_mle[N1] + sigma[1,1]*( p[1] + L_gsg* sigma[1,2] ) + sigma[1,2] *( p[2]+ L_gss * sigma[1,2] )+
0.5 * ( sigma[1,1]**2 * L_ggg + sigma[1,1] *( sigma[1,2]*L_sgg + sigma[2,2]*L_ssg ) + sigma[1,2]*L_sss*sigma[2,2]);
*******************************************************************************;
***************************** Linex Loss Function ****************************;
*******************************************************************************;
U=exp(-z*gamma_mle[N1]); U_g = - z*exp(-z*gamma_mle[N1]); U_gg = z**2 *exp(-z*gamma_mle[N1]);
Exp_Lin_g = U + 0.5* U_gg* sigma[1,1] -
U_g *( sigma[1,1]*( p[1] + L_gsg* sigma[1,2] )+ sigma[1,2]*( p[2] + L_gss* sigma[1,2] ) )-
U_g/2 *( sigma[1,1]**2 *L_ggg+ sigma[1,1]*( sigma[1,2]*L_ssg+ sigma[2,2]*L_ssg)+ sigma[1,2]*L_sss*sigma[2,2] );
*print Exp_Lin_g;
gamma_Lin[N1]= - log((Exp_Lin_g))/z;
*******************************************************************************;
***************************** Precautionary Loss Function ****************************;
*******************************************************************************;
U= gamma_mle[N1]**2; U_g = 2*gamma_mle[N1]; U_gg = 2;
square_pre_g = U + 0.5* U_gg* sigma[1,1] -
U_g *( sigma[1,1]*( p[1] + L_gsg* sigma[1,2] )+ sigma[1,2]*( p[2] + L_gss* sigma[1,2] ) )-
U_g/2 *( sigma[1,1]**2 *L_ggg+ sigma[1,1]*( sigma[1,2]*L_ssg+ sigma[2,2]*L_ssg)+ sigma[1,2]*L_sss*sigma[2,2] );
gamma_Pre[N1]= sqrt((square_pre_g));
/***************************************************************************************;
************************** General Entropy Loss Function *****************************;
***************************************************************************************;
U1_g=-z*gamma_mle[N1]**(-z-1); U1_gg= z*(z+1)*gamma_mle[N1]**(-z-2);
Exp_GE_g = gamma_mle[N1]**(-z)+ 0.5*( U1_gg * gamma_mle[N1]**(-z-2)*sigma[1,1] )-
U1_g *( sigma[1,1]*( p[1]+L_gsg*sigma[1,2] ) +sigma[1,2]*( p[2]+L_gss*sigma[1,2]) )-
U1_g/2 *( sigma[1,1]**2* L_ggg+ sigma[1,1] *( sigma[1,2] *L_sgg+ sigma[2,2]*L_ssg )+sigma[1,2]*L_sss*sigma[2,2] );
gamma_GE[N1]=(Exp_GE_g)**(-1/z);
*********************************************************************************;*/
end;
/*print gamma;
Moment=Gamma_moment[:];
Mle=gamma_mle[:];
SQR=gamma_Sqr[:];
Lin=gamma_Lin[:];
print Moment Mle SQR Lin GE;*/
******************************************************
***************** END Of Simulations ************************;
********************************************************************************************************;
*Bias ;
******;
Bias_mle=abs(gamma_mle[:]-gamma);
Bias_Moment=abs(Gamma_moment[:]-gamma);
Bias_sqr=abs(gamma_Sqr[:]-gamma);
Bias_Lin=abs(gamma_Lin[:]-gamma);
Bias_pre=abs(gamma_pre[:]-gamma);
*MSE ;
******;
MSE_MLE=(gamma_mle-gamma_mle[:])`*(gamma_mle-gamma_mle[:])/k;
MSE_Moment=(Gamma_moment-Gamma_moment[:])`*(Gamma_moment-Gamma_moment[:])/k;
MSE_Sqr=(gamma_Sqr-gamma_Sqr[:])`*(gamma_Sqr-gamma_Sqr[:])/k;
MSE_Lin=(gamma_Lin-gamma_Lin[:])`*(gamma_Lin-gamma_Lin[:])/k;
MSE_pre=(gamma_pre-gamma_pre[:])`*(gamma_pre-gamma_pre[:])/k;
*Relative efficiency;
********************;
Ef_Moment=MSE_MLE/MSE_Moment;
Ef_Sqr=MSE_MLE/MSE_Sqr;
Ef_Lin=MSE_MLE/MSE_Lin;
Ef_pre=MSE_MLE/MSE_pre;
*Relative Bias;
**************;
Rb_MLE=Bias_MLE/gamma*100;
Rb_Moment=Bias_moment/gamma*100;
Rb_Sqr=Bias_Sqr/gamma*100;
Rb_Lin=Bias_Lin/gamma*100;
Rb_pre=Bias_pre/gamma*100;
******************************************;
*PN ;
******;
count1=0; count2=0;count3=0;count4=0;count5=0;count6=0; count7=0; count8=0; count9=0;count10=0;
g_mle=abs(gamma_mle-gamma);
g_Moment=abs(Gamma_moment-gamma);
g_sqr=abs(gamma_Sqr-gamma);
g_Lin=abs(gamma_Lin-gamma);
g_pre= abs(gamma_pre-gamma);
do i=1 to k;
if (g_mle[i] < g_Moment[i]) then count1=count1+1;
if (g_mle[i] < g_sqr[i]) then count2=count2+1;
if (g_mle[i] < g_Lin[i]) then count3=count3+1;
if (g_mle[i] < g_pre[i]) then count7=count4+1;
if (g_Moment[i] < g_sqr[i]) then count4=count5+1;
if (g_Moment[i] < g_Lin[i]) then count5=count6+1;
if (g_Moment[i] < g_pre[i]) then count8=count7+1;
if (g_sqr[i] < g_Lin[i]) then count6=count8+1;
if (g_sqr[i] < g_pre[i]) then count9=count9+1;
if (g_Lin[i] < g_pre[i]) then count10=count10+1;
end;
MLE_VS_Moment= count1/K; MLE_VS_Square= count2/K; MLE_VS_Linex= count3/K; MLE_VS_pre= count4/K;
Moment_VS_Square= count5/K; Moment_VS_Linex= count6/K; Moment_VS_pre= count7/K;
Square_VS_Linex= count8/K; Square_VS_pre= count9/K;
Linex_VS_pre= count10/K;
PN=J(10,1);
PN[1]=MLE_VS_Moment;
PN[2]=MLE_VS_Square;
PN[3]=MLE_VS_Linex;
PN[7]=MLE_VS_pre;
PN[4]=Moment_VS_Square;
PN[5]=Moment_VS_Linex;
PN[8]=Moment_VS_pre;
PN[6]=Square_VS_Linex ;
PN[9]=square_VS_pre;
PN[10]=linex_VS_pre;
print " ";
************************************************************************************;
**** Output Formats for Beta *******;
************************************************************************************;
Bias=J(5,1,0); MSE=J(5,1,0); EF=J(4,1,0);
bias[1]=Bias_mle; bias[2]=Bias_Moment; bias[3]=Bias_sqr; bias[4]=Bias_Lin; bias[5]=Bias_pre;
mse[1]=MSE_MLE; mse[2]=MSE_Moment; mse[3]=MSE_Sqr; mse[4]=MSE_Lin; mse[5]=MSE_pre;
EF[1]=Ef_Moment; EF[2]=Ef_Sqr; EF[3]=Ef_Lin; EF[4]=Ef_pre;
names={MLE, Moment, Square, Linex, pre};
names1={MLE_VS_Moment, MLE_VS_Square, MLE_VS_Linex,MLE_VS_pre, Moment_VS_Square, Moment_VS_Linex,
Moment_VS_pre,Square_VS_Linex, Square_VS_pre, Linex_VS_pre};
print m " " z " " mu1" " mu2" " theta" " alpha1" " alpha2 " " gamma;
print bias [rowname=names] " " mse [rowname=names] " " ;
Print "PN for MLE is better than PN of Moment if PN>0.5" ;
Print PN [rowname=names1] ;
quit;
... View more