Subroutine error

[Salah]

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;

[PaigeMiller]

Show us the entire LOG for this PROC IML so we can see where the error appeared. Please copy the log as text and paste it into the window that appears when you click on the </> icon.

[Rick_SAS]

Looking at the log is essential to understanding errors. In this case, the log tells you 

WARNING: Invalid argument resulted in missing value result.

 count     : number of occurrences is 15
 operation : ## at line 72 column 1
 operands  : y, _TEM1002
Y     15 rows      1 col     (numeric)

_TEM1002      1 row       1 col     (numeric)

 -0.333333

 statement : ASSIGN at line 72 column 1
 traceback : module ADAPTIVE at line 72 column 1

and Line 72 is the statement

72 X = - y ##(-1/del);

The _TEM1002 variable is the temporary expression in parentheses, so apparently
(-1/del) = -0.3333
Apparently, 15 values in the Y vector are nonpositive. The expression 
y##(-1/del)
is undefined when y < = 0 and the power is a floating point number. 

If DEL is an integer, you can rewrite the expression as
X = -1/ x##del;

which should be valid for y ^= 0.

[Salah]

Thank you all for your reply, "del" is not an integer (It could but not necessary).

I am thinking of adding an if statement "If Y<=0 then Y=0.0001" OR "If y<=0 then Go to Step 1" were "Step 1" is defined to take me back to  the beginning of the simulation. However, I tried to change the seed and everything is working fine.  Not sure if I am fooling myself or not. Please advise if any of these ideas are legitimate.

Thank you 

[Rick_SAS · Accepted solution]

> I am thinking of adding an if statement "If Y<=0 then Y=0.0001" OR "If y<=0 then Go to Step 1" 

Because Y is a vector, some elements can be negative whereas others can be positive. Therefore, you must be careful. An IF-THEN statement checks whether ALL elements of vector satisfy the logical condition. See IF-THEN logic with matrix expressions - The DO Loop (sas.com)

It is always better to understand WHY a situation is occurring rather than try to patch it up after it has occurred. I don't know what distribution you are trying to simulate, but I can tell you that the express

y = log ( 1 - (1- U)##(1/theta) ) ;

is negative when U is in the interval (0,1) and theta > 0.

[Salah]

Thank you for your suggestions and help.

I believe I know the reason behind it.  It depends on what data is simulated some of the values of Y are extremely small, not many but even one value can cause a damage!

I noticed that for some seed values I may get no outliers