graph to pick up the parameters for my simulation?

[Salah]

Hi Rick,

I am trying to write a code to simulate data from Bivariate Lomax function with parameters (alpha1, alpha2, C). I will use the created data to compare the different methods of estimation (i.e. MLE, Moment, Bayes estimates such as Square error loss, Linex, Precautionary under different sample sizes and different sets of parameters) based on certain criteria such as Bias, MSE and efficiency.

My concerns are summarized in the followings:

(1) What are the best choices for my parameters? For MLE, Moment and squared error loss I have no problem using almost any values but when it comes to the Linex or Prequationary loss functions, the situation is different....I do get negative values inside the square root or negative values as an argument for the log.

I tried to evaluate the log-likelihood function but it doesn't seem like there is a problem at least not that I am aware of!!

I thought about graphing but I didn't know how and whom to include in the graph.

I would highly appreciate any suggestions.

P.S. The first part (Top part) is the one for evaluating the log-likelihood. My main program is what is shown in green.

Thank you,

Mazen

Proc iml;

seed1=1234;* some times seed=2,0;
seed2=5678;
seed=0;

m=10;


theta=0.1; C=1/theta;


a1=2; a2=3;


********************** Initial values for alpha1 and alpha2 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=prior_g[:];


***********************************************************************;
/*    write the log-likelihood function for bivariate Lomax dist      */
***********************************************************************;

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;

**************************************************************;

X=J(m,1,0);Y=J(m,1,0);

*****************************************************************************;
** Simulating data from bivariate Lomax with parameters alpha1, alpha2 and C ;
*****************************************************************************;

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; ** Simulating X from f(x);
	
	Y[i]=(1+alpha1*X[i])/alpha2*( (1-U2[i])**(-1/(1+c)) -1   );  ** Simulating y from f(Y|x);

end;


/* optional: test the module */

params = {	

		.1 .2 10,
        .1 .5 10,
        .1 0.9 10,

		0.2 0.2 10,
		0.2 0.1 10, 
		
		
		0.5 .1 10,
		0.5 0.2 10,

		0.9 0.1 10,
		0.9 0.2 10,
		0.9 0.5 10,

		1 9 10,
        1 1 10,
        1 0.1 10,
		1 9 40,
        1 1 40,
        	    		  
		1.5 6 10,
        1.5 2 10,
		1.5 1.833 10,
        1.5 9 10,

   		2 1 10,
        2 2 10,
        2 3 10,

		5 45 10,
        5 5 10,
        5 0.56 10,

		6 7.33 10,
        6 1 10,
        6 .667 10,

		10 1 10,
		10 5 10,
		10 8 10
};

LogLik = j(nrow(params),1);
do i = 1 to nrow(params);
   p = params[i,];
   LogLik[i] = LogLik(p);
   
end;



call sort(LogLik); /* ascending by col 1; descending by col 3 */


print Params[c={"alpha1" "alpha2" "C"} label="   "]  "  " LogLik;


quit;

**************** Next is my Main Code ************;

/*Proc iml;

seed1=12345;* some times seed=2,0;
seed2=6301;
seed=0;

m=10;K=1000;


theta=0.1; 


mu1=1.5;mu2=2;


a1=mu1; a2=mu2;

z= 1; 


mu=3;lambda=2 ;  *Initial values for S;

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;


**************************************************************************************;

***************  Optimization Step to find MLEs *********************;

************* Constrain***********************;

    con={1.e-6 1.e-6 1.e-6, . . .};
    opt={1};  
    *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;

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);


*********************************************************************************;


end; * 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;*/

[Rick_SAS]

There seem to be several programs here. I looked at the first one.  It looks like you are sorting the LogLik values, but not the associated parameter values. If you want to know which parameter values correspond to the maximal values of LogLik, try sorting all rows, like this:

results = LogLik || Params;
call sort(results, 1, 1);  /* decreasing by loglik */
print (results[1:8, ])[c={"LogLik" "alpha1" "alpha2" "C"} L="results"]; /* print top 8 */

call scatter(1:nrow(results), results[,1]);

Since your parameter space is three-dimensional, you can't easily graph the LogLiklihood function. However, an attempt is made in the section "Visualizing a function of three variables" in the article "How to find an initial guess for an optimization."  Drawing the picture is not necessary. It is sufficient to sort and print the top parameters. 

[Salah]

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;

[Rick_SAS · Accepted solution]

1. Can you simplify your program so that it shows ONLY the portion that you have a question about? 

2. You might want to look at the article "Trap and cap: Avoid division-by-zero and domain errors when evaluating functions" to see if you can use some of those techniques to avoid domain errors. In an optimization, there are two choices for handling domain errors:

    A. Define a constraint region on which the objective function is well defined.

    B. Make the objective function return a very negative number such as -1e6 if the value of x is not in the domain.

3. Review the "Ten tips before you run an optimization."

[Salah]

Rick you are AWESOMEEEEE thank you so much.

It  works amazingly does!!!!!!!!!!!

I have a slight small question regarding the optimization. In the log page I can still see the error/warning message  (see attachment). I remember that once you said that this shouldn't really concern us if it doesn't happen so often is that true?

Thank you again for your TREMENDOUS help and have a great day

* problem is with either Linex or Precautionary that depends on what initials I am using;


Proc iml;

seed1=12345;* some times seed=2,0;
seed2=6301;
seed=0;

m=10;K=1000;


theta=0.1; 


mu1=1.5;mu2=3;


a1=mu1; a2=mu2;

z= 1; 


mu=3;lambda=2 ;  *Initial values for S;

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;

**************************************************************;
/* Trap bad domain values for the Linex; return missing values when x outside domain */

start Func1(x);
   f = j(nrow(x), ncol(x), .);    /* initialize return values to missing */
   /* trap bad domain values */

   d1 = x;

   domainIdx = loc( d1>0 );     /* find values inside domain */
   if ncol(domainIdx)=0 then return(f); /* all inputs invalid; return missing */
 
   t = x[domainIdx];              /* elements of t are valid */
   numer = log(t);

   f[domainIdx] = numer;  /* evaluate f(t) */
   return( f );
finish;

**************************************************************;
/* Trap bad domain values for the precautionary; return missing values when x outside domain */

start Func2(x);
   f = j(nrow(x), ncol(x), .);    /* initialize return values to missing */
   /* trap bad domain values */

   d1 = x;

   domainIdx = loc( d1>0 );     /* find values inside domain */
   if ncol(domainIdx)=0 then return(f); /* all inputs invalid; return missing */
 
   t = x[domainIdx];              /* elements of t are valid */
   numer = sqrt(t);

   f[domainIdx] = numer;  /* evaluate f(t) */
   return( f );
finish;

******************************************************************;

C_MLE=J(k,1,0); gamma_mle=J(k,1,0);S_mle=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;


**************************************************************************************;

***************  Optimization Step to find MLEs *********************;

************* Constrain***********************;

    con={1.e-6 1.e-6 1.e-6, . . .};
    opt={1};  
    *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];


Gamma_mle[N1]=alpha1_mle/(alpha1_mle+alpha2_mle);
S_mle[N1]=(alpha1_mle+alpha2_mle);




**********************************************************************;


************************************************************************************;
*************     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;

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); 


*******************************************************************************;
*****************************  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);
gamma_Lin[N1]= - Func1(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);

gamma_Pre[N1]= Func2(square_pre_g);



*********************************************************************************;


end; * End Of Simulations;

* problem is with either Linex or Precautionary that depends on what initials I am using;


Linex=gamma_Lin[:];
Precautionary=gamma_Pre[:];
 

print  gamma   Linex "  " Precautionary;




  


quit;

[Rick_SAS]

I'm glad it is working. Hurray!

My advice might or might not apply in your particular circumstance. You need to consider that I have studied your problem for only about 10 minutes. Caveat emptor!

I think what I said is that when you simulate SMALL data samples, the MLE might not converge for a few samples because of random variation that causes the simulated sample to not fit the model.  For LARGE samples, you should not encounter this issue in practice. I have written about how to monitor convergence when you use SAS procedures.  You can do something similar in SAS/IML because the NLPNRA routine returns a code ('rc') which tells you rather the algorithm converged.

In short, I don't know why you sometimes encounter nonconvergence, but by using the 'rc' variable, you can investigate how often it happens and for which parameter values. With some detective work, you might be able to determine the cause. Good luck.

[Salah]

GREAT !!!! Thank you again for everything