BookmarkSubscribeRSS Feed
☑ This topic is solved. Need further help from the community? Please sign in and ask a new question.
Salah
Quartz | Level 8

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;







 

1 ACCEPTED SOLUTION

Accepted Solutions
Rick_SAS
SAS Super FREQ

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.

View solution in original post

5 REPLIES 5
PaigeMiller
Diamond | Level 26

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.

PaigeMiller_0-1663012019648.png

--
Paige Miller
Rick_SAS
SAS Super FREQ

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
Quartz | Level 8

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
SAS Super FREQ

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
Quartz | Level 8

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 

 

 

 

sas-innovate-2024.png

Don't miss out on SAS Innovate - Register now for the FREE Livestream!

Can't make it to Vegas? No problem! Watch our general sessions LIVE or on-demand starting April 17th. Hear from SAS execs, best-selling author Adam Grant, Hot Ones host Sean Evans, top tech journalist Kara Swisher, AI expert Cassie Kozyrkov, and the mind-blowing dance crew iLuminate! Plus, get access to over 20 breakout sessions.

 

Register now!

Multiple Linear Regression in SAS

Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.

Find more tutorials on the SAS Users YouTube channel.

From The DO Loop
Want more? Visit our blog for more articles like these.
Discussion stats
  • 5 replies
  • 720 views
  • 0 likes
  • 3 in conversation