Hello I am using my previous code which is working fine, just changed the density. I keep getting warning messages "WARNING: Invalid argument resulted in missing value result." I check the line that caused this issue and printed it to see all values, but I found no problem with it. Can you please help me figure out it and solve this issue. Thank you proc iml;
*seed=0; NN=10;
theta1=2; theta2=10; del=3;
case=2; n1=30; m1=15; k=J(m1,1,0); K[m1]=n1-m1; L=K; n2=n1; m2=m1;
**********************************************************************************************************************;
*********************** Exact OVL *******************************;
R=(theta1/theta2);
R_star=((m2-1)/m2)*R;
R1=R_star;
rho_Exact=2*sqrt(R1)/(1+R1);
lambda_Exact=(2*sqrt(R1)/(1+R1))**2;
delta_Exact=1-R1**(1/(1-R1))*abs(1-1/R1);
******************************************************************************;
start Uniform_p(m,R) global(seed);
V=J(m,1,0); U=J(m,1,0);W=J(m,1);
W=randfun(m,"Uniform");
Do i=1 to m;
Sum=0;
Do j=m-i+1 to m;
Sum=Sum+R[j];
End;
Sum=Sum+i;
V[i]=W[i]**(1/Sum);
End;
Do i=1 to m;
U[i]=1-prod ( V[m: (m-i+1)] );
End;
Call sort (U);
return U;
Finish;
********************************************************************;
Start Adaptive(n,m,R,theta) Global(del) ;
U=uniform_p(m,R);
Y=( log ( 1 - (1- U)##(1/theta) ) ) ;
X = - y ##(-1/del);
T=X[m]+2;
Do idx=1 to m-1;
If (( X[idx] < T) & (T <= X[idx+1] )) then j=idx;
End;
If X[m]>T then
Do;
W=Uniform_p(m-j,R);
WE= (1 - Exp(-1/T##del));
X_m_j =( -1/ ( log ( 1 - ( (1- W)##(1/theta) ) # WE ) ) )##del;
X_Adptive= X[1:j]//X_m_j;
Rm=n-m-sum(R[1:J]);
R2=J(m-j,1,0);R2[m-J]=Rm;
newR=R[1:j]//R2;
End;
Else
Do;
j= m;
newR=R;
X_Adptive=X;
End;
P1=[X_Adptive,j,newR,y];
return P1;
Finish;
********** Subroutine to compute the MLEs **************;
*********************************************************;
start MLE(X,n,m,R,J);
Part=1-Exp(- 1/X##(-del));
deno = R[m]*log(part[m])+ sum (log(part)) + sum( (R#log(part))[1:J] ) ;
Theta_MLE = (- m/deno);
return(Theta_MLE);
finish;
******************************************************************************************;
theta1_Adp=J(NN,1,0); theta2_Adp=J(NN,1,0);H_Adp=J(NN,1,0);
R_star_Adp=J(NN,1,0);
rho_star_Adp=J(NN,1,0);Var_rho_Adp=J(NN,1,0);Bias_rho_Adp=J(NN,1,0);
lambda_star_Adp=J(NN,1,0);Var_lambda_Adp=J(NN,1,0);Bias_lambda_Adp=J(NN,1,0);
delta_star_Adp=J(NN,1,0);Var_delta_Adp=J(NN,1,0);Bias_delta_Adp=J(NN,1,0);
MSE_rho_Adp=J(NN,1,0);MSE_lambda_Adp=J(NN,1,0);MSE_delta_Adp=J(NN,1,0);
Lower_rho_Adp=J(NN,1);Upper_rho_Adp=J(NN,1);Length_rho_Adp=J(NN,1);
Lower_lambda_Adp=J(NN,1);Upper_lambda_Adp=J(NN,1);Length_lambda_Adp=J(NN,1);
Lower_delta_Adp=J(NN,1);Upper_delta_Adp=J(NN,1);Length_delta_Adp=J(NN,1);
COV1_Adp=J(NN,1,0);COV2_Adp=J(NN,1,0);COV3_Adp=J(NN,1,0);***COVi= stands for coverage of OVLi, i=1,2,3**;
var_rho_Adp=J(NN,1,0);var_lambda_Adp=J(NN,1,0);var_delta_Adp=J(NN,1,0);
MSE1=J(NN,1,0);MSE2=J(NN,1,0);MSE3=J(NN,1,0);
Do JJ=1 to NN;
X=J(m1,1,0);Y=J(m2,1,0); X_ADP=J(m1,1,0);Y_ADP=J(m2,1,0);
*********************************************************************;
AdaptiveResX=Adaptive(n1,m1,K,theta1);
X_ADP=AdaptiveResX$1;
J1=AdaptiveResX$2;
newK=AdaptiveResX$3;
y=AdaptiveResX$4;print y;
K=newK;
AdaptiveResY=Adaptive(n2,m2,L,theta2);
Y_ADP=AdaptiveResY$1;
J2=AdaptiveResY$2;
newL=AdaptiveResY$3;
L=newL;
Theta1_Adp[JJ] = MLE(X_ADP,n1,m1,K,J1);
Theta2_Adp[JJ] = MLE(Y_ADP,n2,m2,L,J2);
R_star_Adp[JJ]=(m2-1)/m2*(theta1_Adp[JJ]/theta2_Adp[JJ]);
End;
** The estimated values of the OVL **;
rho_star_Adp=2#sqrt(R_star_Adp)#(1+R_star_Adp)##(-1);
lambda_star_Adp=rho_star_Adp##2;
delta_star_Adp=1-R_star_Adp##(1/(1-R_star_Adp))#abs(1-1/R_star_Adp);
****** Asymptotic Variance of the Overlap ********;
var_rho_Adp=R_star_Adp#(1-R_star_Adp)##2/(1+R_star_Adp)##4#(m1+m2-1)/(m1*(m2-2));
var_lambda_Adp=16#R_star_Adp##2#(1-R_star_Adp)##2/(1+R_star_Adp)##6#(m1+m2-1)/(m1*(m2-2));
var_delta_Adp=(m1+m2-1)/(m1*(m2-2))#R_star_Adp##(2/(1-R_star_Adp))#(log(R_star_Adp))##2/(1-R_star_Adp)##2;
****** Asymptotic Bias of the Overlap ********;
Bias_rho_Adp=(m1+m2-1)/( m1*(m2-2))# sqrt(R_star_Adp)# ( 3#R_star_Adp##2-6#R_star_Adp-1)/( 4#(1+R_star_Adp)##3 );
Bias_lambda_Adp=(m1+m2-1)/(m1*(m2-2))#4#R_star_Adp##2#(R_star_Adp-2)/(1+R_star_Adp)##4;
H_Adp= R_star_Adp##(1/(1-R_star_Adp)) # ( R_star_Adp#(2#R_star_Adp-log(R_star_Adp)-2)#log(R_star_Adp) - (R_star_Adp-1)##2 )/(R_star_Adp-1)##3;
H_Adp=R_star_Adp##2 # R_star_Adp##((2#R_star_Adp-1)/(1-R_star_Adp)) # ( R_star_adp #
(2#R_star_Adp-log(R_star_Adp)-2)#log(R_star_Adp) - (R_star_Adp-1)##2 )/(R_star_Adp-1)##3;
Do JJ=1 to NN;
if (R_star_Adp[JJ] > 1) then Bias_delta_Adp[JJ]=(m1+m2-1)/(2*m1*(m2-2))* H_Adp[JJ]/2;
if (R_star_Adp[JJ] < 1) then Bias_delta_Adp[JJ]=-(m1+m2-1)/(2*m1*(m2-2))* H_Adp[JJ]/2;
End;
Theta1_h=theta1_Adp[:];
Theta2_h=theta2_Adp[:];
Bias_rho=Bias_rho_Adp[:]; Bias_lambda=Bias_lambda_Adp[:]; Bias_delta=Bias_delta_Adp[:];
** I used MSE below b/c it provides smaller MSE values;
MSE1 = ( rho_star_Adp - rho_Exact )##2;
MSE2 = ( lambda_star_Adp - lambda_exact)##2;
MSE3 = ( delta_star_Adp - delta_exact)##2;
Var1= Var_rho_Adp[:];
Var2= Var_lambda_Adp[:];
Var3= Var_delta_Adp[:];
**** Interval estimation using Asymptotic technique ****;
Do JJ=1 to NN;
*************************** OVL1 *****************************;
Lower_rho_Adp[JJ] = rho_star_Adp[JJ] - Bias_rho_Adp[JJ] - 1.96*sqrt(Var_rho_Adp[JJ]);
Upper_rho_Adp[JJ] = rho_star_Adp[JJ] - Bias_rho_Adp[JJ] + 1.96*sqrt(Var_rho_Adp[JJ]);
Length_rho_Adp[JJ]=Upper_rho_Adp[JJ]-Lower_rho_Adp[JJ];
If ( (rho_exact >= Lower_rho_Adp[JJ])& (rho_exact <= upper_rho_Adp[JJ]) ) then COV1_Adp[JJ]=1;* COV is the coverage;
*************************************************************************************;
*************************** OVL2 *****************************;
Lower_lambda_Adp[JJ] = lambda_star_Adp[JJ] - Bias_lambda_Adp[JJ] - 1.96*sqrt(Var_lambda_Adp[JJ]);
Upper_lambda_Adp[JJ] = lambda_star_Adp[JJ] - Bias_lambda_Adp[JJ] + 1.96*sqrt(Var_lambda_Adp[JJ]);
Length_lambda_Adp[JJ]= Upper_lambda_Adp[JJ] - Lower_lambda_Adp[JJ];
If ( (lambda_exact >= Lower_lambda_Adp[JJ])& (lambda_exact <= upper_lambda_Adp[JJ]) )
then COV2_Adp[JJ]=1;
**************************************************************************************;
*************************** OVL3 ******************************;
Lower_delta_Adp[JJ] = delta_star_Adp[JJ] - Bias_delta_Adp[JJ] - 1.96*sqrt(Var_delta_Adp[JJ]);
Upper_delta_Adp[JJ] = delta_star_Adp[JJ] - Bias_delta_Adp[JJ] + 1.96*sqrt(Var_delta_Adp[JJ]);
Length_delta_Adp[JJ]= Upper_delta_Adp[JJ] -Lower_delta_Adp[JJ];
If ( (delta_exact >= Lower_delta_Adp[JJ])& (delta_exact <= upper_delta_Adp[JJ]) )
then COV3_Adp[JJ]=1;
end;
***************** END Of Simulations ************************;
********************************************************************************************************;
BiasOVL1_Adp=ABS(Bias_rho_Adp[:]); ** why I added 0.001 is explained above ;
BiasOVL2_Adp=ABS(Bias_lambda_Adp[:]);
BiasOVL3_Adp=ABS(Bias_delta_Adp[:]);
MSE11=MSE1[:];
MSE22=MSE2[:];
MSE33=MSE3[:];
length_OVL1_Adp=Length_rho_Adp[:];
length_OVL2_Adp=Length_lambda_Adp[:];
length_OVL3_Adp=Length_delta_Adp[:];
coverage1_Adp=COV1_Adp[:];
coverage2_Adp=COV2_Adp[:];
coverage3_Adp=COV3_Adp[:];
****formatting the output *******;
*print Lower_rho_RSS exact_rho upper_rho_RSS;
Bias=J(3,1);Length=J(3,1);coverage=J(3,1); MSE=J(3,1,0);
Bias[1,1]=(BiasOVL1_Adp ); Bias[2,1]=(BiasOVL2_Adp ); Bias[3,1]=(BiasOVL3_Adp );
MSE[1,1]= MSE11; MSE[2,1]= (MSE22 ); MSE[3,1]= (MSE33 );
Length[1,1]=(length_OVL1_Adp ); Length[2,1]=(length_OVL2_Adp ); Length[3,1]=(length_OVL3_Adp );
coverage[1,1]=(COV1_Adp[:]); coverage[2,1]=(COV2_Adp[:] ); coverage[3,1]=(COV3_Adp[:] );
names={rho,lambda,delta};
Print "Results from T3";
print case n1 m1;
print names bias " " MSE " " Length " " coverage ;
quit;
... View more
