Subroutine error

Salah — Wed, 13 Sep 2023 19:47:49 GMT

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;

Re: Subroutine error

PaigeMiller — Wed, 13 Sep 2023 19:53:46 GMT

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.

Re: Subroutine error

Rick_SAS — Wed, 13 Sep 2023 20:02:21 GMT

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.

Re: Subroutine error

Salah — Thu, 14 Sep 2023 05:58:36 GMT

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

Re: Subroutine error

Rick_SAS — Thu, 14 Sep 2023 09:50:18 GMT

> 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.

Re: Subroutine error

Salah — Thu, 14 Sep 2023 16:59:08 GMT

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

topic Subroutine error in SAS/IML Software and Matrix Computations

Subroutine error

Re: Subroutine error

Re: Subroutine error

Re: Subroutine error

Re: Subroutine error

Re: Subroutine error