Hello I am trying to use SAS/IML to solve a system of nonlinear equations. In the beginning, I succeeded till I used smaller values for the parameters, then I start receiving error messages (see attached log file). I noticed that the errors only occur when I do more simulations, 1000 for example. No errors when I run 1 or 2 simulations. I also want to say that I read your answer to a similar problem (see the link below) but I couldn't implement it to my code. Didn't know where to put opt? There was no example in the link so I couldn't use that help. https://communities.sas.com/t5/SAS-IML-Software-and-Matrix/How-to-get-the-optimizations-results-when-the-error-occurs/td-p/206133 Thank you %macro ODSOff(); ods graphics off; ods exclude all; ods noresults; %mend; %macro ODSOn(); ods graphics on; ods exclude none; ods results; %mend; %odsoff; proc iml; seed=12345; k=1000; n=16; alpha =0.7; Lambda=1.5; *be careful for PWM, I had to assume that Lambda >1; m=n; R=J(m,1,0); **** MLE ***; start MLE_func(y) global (n,m,X,R); func=0; alpha=y[1]; lambda=y[2]; Sum_log=J(m,1,0); Sum_log=log(x); Sum_log_1=J(m,1,0); Sum_log_1=log(1-X##lambda); func=func + m*log(alpha*lambda)+(lambda-1)* Sum_log[+] + (alpha-1) * Sum_log_1[+] ; Return(func); finish; *** Solve for MOM: by solving system of equations ***; *******************************************************; /* I solved the system of nonlinear equations using: https://blogs.sas.com/content/iml/2018/02/28/solve-system-nonlinear-equations-sas.html */ start MOM_Fun(var) global (n,m,X,R); X_bar=X[:]; Xsqr=X##2; alpha = var[1]; lambda = var[2]; f = j(1, 2, .); /* return a ROW VECTOR */ f[1] = X_bar-gamma(alpha+1)*gamma(1/lambda+1)/gamma(alpha+1/lambda+1); f[2] = Xsqr[:] - alpha*gamma(alpha)*gamma(2/lambda+1)/gamma(alpha+2/lambda+1); return (f); finish; *** Solve for PWM: by solving system of equations ***; *******************************************************; start PWM_Fun(var) global (n,m,X,R); alpha = var[1]; lambda = var[2]; call sort (X); X_bar=X[:]; Sum=0; do i=1 to n; sum=sum+(n-i)*X[i]; end; sum1=sum/(n*(n-1)); pi=constant("pi"); B=-lambda*gamma(1-1/lambda); f = j(1, 2, .); f[1] = X_bar + pi* alpha * CSC(pi/lambda)* gamma(alpha)/( gamma( alpha + 1/lambda + 1)* B ); f[2] = sum1 + pi* alpha * CSC(pi/lambda)* gamma(2*alpha)/( gamma( 2*alpha + 1/lambda + 1)* B ); return (f); finish; *** Solve for SDPWM: by solving system of equations ***; *******************************************************; start SDPWM_Fun(var) global (n,m,X,R); alpha = var[1]; lambda = var[2]; call sort (X); X_bar=X[:]; Sum= (1-X##lambda)##alpha # X; pi=constant("pi"); B=-lambda*gamma(1-1/lambda); f = j(1, 2, .); f[1] = X_bar + pi* alpha * CSC(pi/lambda)* gamma(alpha)/( gamma( alpha + 1/lambda + 1)* B ); f[2] = sum[:] + pi* alpha * CSC(pi/lambda)* gamma(2*alpha)/( gamma( 2*alpha + 1/lambda + 1)* B ); return (f); finish; *** Solve for Shanon Entropy: by solving system of equations ***; *******************************************************; start Shanon_Fun(var) global (n,m,X,R); alpha = var[1]; lambda = var[2]; call sort (X); LNX= log(X); LNXB=Log(1-X##lambda); Sum=LNX#digamma(1-X##lambda); f = j(1, 2, .); f[1] = LNX[:]+( digamma(alpha+1)- digamma(1) )/lambda; f[2] = LNXB[:] + ( digamma(alpha+1)- digamma(alpha) ); *f[3]= (1-alpha)*Sum[:]+( digamma(alpha+1)- digamma(1) -1 )/lambda; return (f); finish; ********************* Simulation Step ***************************; alpha_MLE=J(k,1,0); lambda_MLE=J(k,1,0); alpha_Mom=J(k,1,0); lambda_Mom=J(k,1,0); alpha_PMW=J(K,1,0); lambda_PMW=J(K,1,0); alpha_SDPMW=J(k,1,0); lambda_SDPMW=J(k,1,0); alpha_Shanon=J(k,1,0); lambda_Shanon=J(k,1,0); MSE_MLE_a=J(K,1,0); MSE_MLE_L=J(K,1,0); MSE_MOM_a=J(K,1,0); MSE_MOM_L=J(K,1,0); MSE_PMW_a=J(K,1,0); MSE_PMW_L=J(K,1,0); MSE_SDPMW_a=J(K,1,0); MSE_SDPMW_L=J(K,1,0); MSE_Shanon_a=J(K,1,0); MSE_Shanon_L=J(K,1,0); Do N1=1 to K; W=J(m,1,0);V=J(m,1,0);X_p=J(m,1,0); X=J(m,1,0); U=J(m,1,0); do i=1 to m; W[i]=Uniform(seed); end; S=1+R[m]; V[1]=W[1]**(1/S); do i=2 to m; S=S+(1+R[m-i+1]); V[i]=W[i]**(1/S); end; do i=1 to m; U[i]=1-prod( V[ m:( m-i+1)] );** the U's are the required progressively type II from U(0,1); END; X=( 1-(1-U)##(1/alpha) )##(1/lambda); ****************************************************************** ; *** MLE for the Kumaraswamy ***** ; *******************************************************************; ************* Constrain MLE ***********************; con={0.01 0.01 , . .}; x0={0.05, 0.05}; opt={2 1}; tc={10000 14000}; call nlpnra(rc, MLE_ret, "MLE_func", x0, opt, con,tc); alpha_mle[N1] =MLE_ret[1]; lambda_mle[N1]=MLE_ret[2]; ************************* Find MOM ************************; /* x[1] x[2] constraints. Lower bounds in 1st row; upper bounds in 2nd row */ con = {1e-3 1e-3, . .}; /* x[1] > 0 and x[2] > 0; */ x0 = {0.5 1.1}; /* initial guess */ optn = {2 /* solve least square problem that has 3 components */ 1}; /* amount of printing */ call NLPLM(rc, Soln, "MOM_Fun", x0, optn) blc=con; /* or use nlphqn */ alpha_mom[N1]=Soln[1]; lambda_mom[N1]=Soln[2]; ************************* Find PWM ************************; call nlphqn(rc, Soln1, "PWM_Fun", x0, optn) blc=con; alpha_PMW[N1]=Soln1[1]; lambda_PMW[N1]=Soln1[2]; ************************ Find SDPWM ************************; call nlphqn(rc, Soln2, "SDPWM_Fun", x0, optn) blc=con; alpha_SDPMW[N1]=Soln2[1]; lambda_SDPMW[N1]=Soln2[2]; ************************* Find Shanon ************************; call nlphqn(rc, Soln3, "Shanon_Fun", x0, optn) blc=con; alpha_Shanon[N1]=Soln3[1]; lambda_Shanon[N1]=Soln3[2]; *MSE alpha; **********; MSE_MLE_a[N1]=(alpha_mle[N1]-alpha)**2; MSE_MOM_a[N1]=(alpha_MOM[N1]-alpha)**2; MSE_PMW_a[N1]=(alpha_PMW[N1]-alpha)**2; MSE_SDPMW_a[N1]=(alpha_SDPMW[N1]-alpha)**2; MSE_Shanon_a[N1]=(alpha_Shanon[N1]-alpha)**2; *MSE Lambda; **********; MSE_MLE_L[N1]=(Lambda_mle[N1]-Lambda)**2; MSE_MOM_L[N1]=(Lambda_MOM[N1]-Lambda)**2; MSE_PMW_L[N1]=(Lambda_PMW[N1]-Lambda)**2; MSE_SDPMW_L[N1]=(Lambda_SDPMW[N1]-Lambda)**2; MSE_Shanon_L[N1]=(Lambda_Shanon[N1]-Lambda)**2; **************************************************************************; End; %odson MLE_alpha_hat=alpha_mle[:]; MLE_lambda_hat=lambda_mle[:]; Mom_alpha_hat =alpha_Mom[:]; Mom_lambda_hat=lambda_Mom[:]; PMW_alpha_hat =alpha_PMW[:]; PMW_lambda_hat=lambda_PMW[:]; SDPMW_alpha_hat =alpha_SDPMW[:]; SDPMW_lambda_hat=lambda_SDPMW[:]; Shanon_alpha_hat =alpha_Shanon[:]; Shanon_lambda_hat=lambda_Shanon[:]; print alpha lambda; print MLE_alpha_hat MLE_lambda_hat; print Mom_alpha_hat Mom_lambda_hat; Print PMW_alpha_hat PMW_lambda_hat; Print SDPMW_alpha_hat SDPMW_lambda_hat; Print Shanon_alpha_hat Shanon_lambda_hat; quit;
... View more