Hi Rick, I fixed the Fisher matrix, it is positive definite and symmetric. Unfortunately this didn't help with the values of Exp_Lin_g. I still get negative values. Exp_Lin_g must be a number between 0 and 1 as the estimate must be >0, hence - Log(Exp_Lin_g)>0 Upon checking my code carefully I was surprised to find the following: (1) I get troubles with certain choices of the seed (negative values for Exp_Lin_g) can't tell which seeds as they keep changing!!!! (2) I couldn't explain or understand how choosing the initials that provide the highest value for the likelihood creates errors while choices such as 1.5 and 3 which do not maximize the likelihhod function won't really work ( i.e. Exp_Lin_g<0, Exp_Lin_g depends heavily on the MLE). Proc iml;
seed1=123;seed2=456;seed=0;
m=20; C=10;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;
alpha1=Prior_a1[:]; alpha2=Prior_a2[:]; 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;
*****************************************************************************;
** Simulating data from bivariate Lomax with parameters alpha1, alpha2 and C ;
*****************************************************************************;
X=J(m,1,0);Y=J(m,1,0);
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.5 3 10, 0.9 0.9 10, 0.9 .1 10, 0.1 0.9 10, 0.7 0.7 10, 0.7 1 10,4 6 10,5 4 10, 5.5 4 10,
1 0.7 10, 1.5 5 10, 1.5 5.5 10, 1.5 6 10, 5 1.5 10, 5.5 1.5 10, 6 1.5 10,5 5 10,5 5.5 10,
2 2 10, 2 5 10, 2 5.5 10, 2 6 10, 5 2 10, 5.5 2 10, 6 2 10, 6 4 10, 5 6 10,
3 5 10, 3 5.5 10, 3 6 10, 5 3 10, 5.5 3 10, 6 3 10,4 5 10,4 5.5 10,5.5 5 10,6 5 10};
LogLik = j(nrow(params),1);
do i = 1 to nrow(params);
p = params[i,];
LogLik[i] = LogLik(p);
end;
results = LogLik || Params;
call sort(results, 1, 1); /* decreasing by loglik */
print (results[1:30, ])[c={"LogLik" "alpha1" "alpha2" "C"} L="results"];
call scatter(1:nrow(results), results[,1]);
quit;
(3) Increasing the stimulation size results on increasing the chances of getting awful values for Exp_Lin_g...I mean getting large negative/positive values such as -2.12 or 3.7!!!!! (4) Replacing negative values with missing values is not working. I called it "gamma_lin_Rick and as you see in the pic below it is negative and it supposed to be "."!!!! I am attaching my code below. Thank you Rick. Proc iml;
*seed1=2; *seed2=22; *These choices are working sometimes;
*seed1=123; *seed2=0; seed1=123; seed2=456;
m=10;K=1000;
a1=1.5;a2=2;
a1=1.5; a2=3; /*this choice caused mistakes*/
C=10;z= 1;mu=3;lambda=2 ; *Initial values for S;
* 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]=RAND('GAMMA',a1);
Prior_a2[i]=RAND('GAMMA',a2);
prior_g[i]= RAND('BETA',a1,a2);
end; *DO LOOP OF THE GAMMA;
alpha1=Prior_a1[:]; *true value for alpha1;
alpha2=Prior_a2[:]; *true value for alpha2;
gamma= prior_g[:];
***********************************************************************;
*** Log likelihood function
************************************************************************;
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; return missing values when x outside domain */
start Func1(x);
f = j(nrow(x), ncol(x), .); /* initialize return values to missing */
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;
/* The part below is for simulating data from Bivariate Lomax*/
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_Lin_Rick=J(k,1,0);Exp_Lin_g=J(k,1,0);
X=J(m,1,0);Y=J(m,1,0);
Do N1=1 to K;
U1=J(m,1,0);U2=J(m,1,0);
do i=1 to m;
U1[i]=Uniform(seed1);
U2[i]=Uniform(seed2);
/* the method below is a better way for simualtion,I wanted to control the seed so I can figure out what went wrong!!!!!
/* But I failed!!!!!!
*U1[i]=RAND('UNIFORM');
*U2[i]=RAND('UNIFORM');
*/
X[i]=( (1-U1[i])**(-1/c) -1 )/alpha1;
Y[i]=(1+alpha1*X[i])/alpha2*( (1-U2[i])**(-1/(1+c)) -1 );
end;
/* End of the simulation */
*************** Optimization Step to find MLEs *********************;
************* Constrain***********************;
con={1.e-6 1.e-6 1.e-6, . . .};
opt={1};
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[+]) ;
**************************************************************************;
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=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]= (2*m)/S_MLE[N1]**2 -( 2+c_mle[N1] )* Sum_I22[+];
sigma=inv(fisher);
*print (eigval(sigma)); /*this is to check is sigma is positive definite*/
*******************************************************************************;
***************************** 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]);
x1=0.5* U_gg* sigma[1,1];
x2=sigma[1,1]*( p[1] + L_gsg* sigma[1,2] );
x3=sigma[1,2]*( p[2] + L_gss* sigma[1,2] );
x4=sigma[1,1]**2 *L_ggg;
x5=sigma[1,1]*( sigma[1,2]*L_ggs+ sigma[2,2]*L_ssg);
x6=sigma[1,2]*L_sss*sigma[2,2];
Exp_Lin_g[N1] = U + x1 + U_g *( x2 + x3 )+ 0.5 *U_g *( x4+ x5+ x6 );
gamma_Lin_Rick[N1]= -Func1(Exp_Lin_g[N1])/z;
gamma_Lin[N1]= -Log(Exp_Lin_g[N1])/z;
end; * End Of Simulations;
Estimate_Linex=gamma_Lin[:];
MSE=(gamma_Lin-gamma_Lin[:])`*(gamma_Lin-gamma_Lin[:])/k;
print Estimate_Linex " " MSE;
print gamma_Lin gamma_Lin_Rick ;
quit;
... View more
