Hello My code is supposed to simulate two samples from exponential distribution based on progressive type II. Then, I must modify the samples by truncating each one at a particular time and replacing the truncated values with new ones. The final samples should be used to estimate the overlap function between the two samples. I have three truncation choices: after 30%, 80%, or 0% of the data. My problem is that all results are almost identical. No matter if I do the truncation or not!!! I also found that if I fix the percentage of the truncation and change the censoring scheme (i.e., to censor the data at the beginning, at the end, or in the middle(results remain unchanged) which also doesn't make any sense!!! I would greatly appreciate any help in resolving this matter. proc iml;
*Results for T1, T2 and T3 have striking similarities, I checked everything and don't know what is the problem,
Results for Ti, i=1,2,3 when changing the censoring scheme also produce stricking similarities!!!!
**********************************************************************************************************************;
seed=0; NN=10000;*NN is the number of simulations;
theta1=1; theta2=5;
case=17; n1=80; m1=50; k=J(m1,1,0); K[m1]=n1-m1; L=K; n2=n1; m2=m1;
/**********************************************************************************************************************;
case=16; n1=80; m1=50; k=J(m1,1,0); K[1]=n1-m1; L=K; n2=n1; m2=m1;
case=17; n1=80; m1=50; k=J(m1,1,0); K[m1]=n1-m1; L=K; n2=n1; m2=m1;
case=18; n1=80; m1=50; k=J(m1,1,0); k[41:50]=3; L=K; n2=n1; m2=m1;
*/
*********************** Exact OVL *******************************;
R=(theta1/theta2);
R_star=((m2-1)/m2)*R;
Exact_rho=2*sqrt(R_star)/(1+R_star);
Exact_lambda=(Exact_rho)**2;
Exact_delta=1-R_star**(1/(1-R_star))*abs(1-1/R_star);
******************************************************************************;
start Uniform_p(m,R) global(seed);
** The Required Progressively Type-II from U(0,1);
V=J(m,1,0); U=J(m,1,0);W=J(m,1);
Call randseed(seed);
Call randgen(W, "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(X, T, n,m,R,theta) ;
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);
X_m_j = T - theta*log(1-W);
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];
return P1;
Finish;
********** Subroutine to compute the MLEs **************;
*********************************************************;
start MLE(X,n,m,R,Rm,J);
Theta_MLE=( Rm*X[m]+X[+]+ sum( (R#X)[1:J] ) )/m;
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);
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);
**** Data set I;
***************;
X = - theta1 * log (1- uniform_p(m1,K)); * Create Progressive type II data;
*T1=X[ceil(0.3*m1)]; * Keep the top 80% and replace the remaining data using Adaptive module;
*T1=X[ceil(0.8*m1)]; * Keep the top 30% and replace the remaining data using Adaptive module;
T1=X[m1]+2; * Do not truncate the data;
AdaptiveResX=Adaptive(X, T1,n1,m1,K,theta1);**Create Adaptive hybrid progressive type II;
X_ADP=AdaptiveResX$1;* This is my Adaptive IIPH that I will use to find MLEs;
J1=AdaptiveResX$2;
newK=AdaptiveResX$3;
K=newK; Km=K[m1];
**** Data set II;
****************;
Y = - theta2 * log (1- uniform_p(m2,L));
*T2=Y[ceil(0.3*m2)];
*T2=Y[ceil(0.8*m2)];
T2=Y[m2]+2;
AdaptiveResY=Adaptive(Y, T2,n2,m2,L,theta2);
Y_ADP=AdaptiveResY$1;* This is my Adaptive IIPH that I will use to find MLEs;
J2=AdaptiveResY$2;
newL=AdaptiveResY$3;
L=newL; Lm=L[m2];
*print J1 T1 J2 T2 k newk l newl;
*print X T1 J1 X_adp " " Y T2 J2 Y_adp;
*********************************************************************;
*Theta1_Adp[JJ]=( Km* X_ADP[m1]+X_ADP[+] + sum( ( K#X_ADP )[1:J1] ) )/m1;
*Theta2_Adp[JJ]=( Lm* Y_ADP[m2]+Y_ADP[+] + sum( ( L#Y_ADP )[1:J2] ) )/m2;
Theta1_Adp[JJ]=MLE(X_ADP,n1,m1,K,Km,J1);
Theta2_Adp[JJ] =MLE(Y_ADP,n2,m2,L,Lm,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 Overlape ********;
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=R_star_Adp##(2/(1-R_star_Adp))#(log(R_star_Adp))##2/(1-R_star_Adp)##2#(m1+m2-1)/(m1*(m2-2));
****** Asymptotic Bias of the Overlape ********;
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;
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];
if (R_star_Adp[JJ] < 1) then Bias_delta_Adp[JJ]=-(m1+m2-1)/(2*m1*(m2-2))* H_Adp[JJ];
end;
Theta1_h=theta1_Adp[:];
Theta2_h=theta2_Adp[:];
Bias_rho=Bias_rho_Adp[:]; Bias_lambda=Bias_lambda_Adp[:]; Bias_delta=Bias_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 ( (Exact_rho >= Lower_rho_Adp[JJ])& (Exact_rho <= 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 ( (Exact_lambda >= Lower_lambda_Adp[JJ])& (Exact_lambda <= 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 ( (Exact_delta >= Lower_delta_Adp[JJ])& (Exact_delta <= upper_delta_Adp[JJ]) )
then COV3_Adp[JJ]=1;
end;
***************** END Of Simulations ************************;
********************************************************************************************************;
BiasOVL1_Adp=ABS(Bias_rho_Adp[:]);
BiasOVL2_Adp=ABS(Bias_lambda_Adp[:]);
BiasOVL3_Adp=ABS(Bias_delta_Adp[:]);
varOVL1=VAR_rho_Adp[:];
VAROVL2=VAR_lambda_Adp[:];
VAROVL3=VAR_delta_Adp[:];
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[:];
****formating the output *******;
*print Lower_rho_RSS exact_rho upper_rho_RSS;
Bias=J(3,1);Length=J(3,1);coverage=J(3,1); Var=J(3,1,0);
bias[1,1]=Round(BiasOVL1_Adp,0.0001 ); bias[2,1]=Round(BiasOVL2_Adp,0.0001 ); bias[3,1]=Round(BiasOVL3_Adp,0.0001 );
Var[1,1]=Round(VarOVL1,0.0001 ); Var[2,1]=Round(VarOVL2,0.0001 ); Var[3,1]=Round(VarOVL3,0.0001 );
Length[1,1]=Round(length_OVL1_Adp,0.0001 ); Length[2,1]=Round(length_OVL2_Adp,0.0001 ); Length[3,1]=Round(length_OVL3_Adp,0.0001 );
coverage[1,1]=Round(COV1_Adp[:],0.001 ); coverage[2,1]=Round(COV2_Adp[:],0.001 ); coverage[3,1]=Round(COV3_Adp[:],0.001 );
names={rho,lambda,delta};
print case R n1 m1;
print names bias " " Var " " Length " " coverage ;
quit;
... View more
