LEVMAR Optimization cannot be completed!!!

[Salah]

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

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;

[Rick_SAS]

There is a general way to solve these problems.

1. Scroll down in the log until you see the first ERROR message. In this case, it says

ERROR: (execution) Invalid argument to function.operation :
GAMMA at line 76 column 37
operands : _TEM1004_TEM1004 1 row 1 col (numeric)1000001

This tells you that the problem is in a call to the GAMMA function, which is on Line 76 (Col=37) of the log. It says that the operand (which means argument to the function) is a 1x1 matrix that has the value 1000001.

2. Go to line 76 in the log. There you will see the call to the GAMMA function. I don't know which call is failing, but it might be GAMMA(1/lambda+1). Notice that lambda is not constrained away from zero, so it is possible that you are seeing lambda->0, which causes the GAMMA function of 1/lambda to blow up.

So, what to do?  Your program is too long (275 lines) for me to try to understand it, so my suggestion is that you try to isolate the issue and reproduce it in a simpler program. When you isolate the problem, the solution will sometimes be more apparent. Some other ideas:

1. In optimization problems, it is often sufficient to optimize the LOG of a function. Can you work with the LOG of the expression? LOG( -gamma(alpha+1)*gamma(1/lambda+1)/gamma(alpha+1/lambda+1)) can be expressed in terms of the 

LGAMMA function.

2. Be sure you can solve ONE problem before you attempt to embed the problem in a simulation study that needs to solve thousands of problems. 

3. Explain in words why you are trying to solve these simultaneous equations. There might be an easier way. 

[Rick_SAS]

Also, please double-check your second equation in the MOM function. You seem to be using the second (non-centered) moment. Did you intend to use the variance or the second centered moment? I've only used this method to match the mean and the variance.

[Salah]

Thank you for the reply. I am trying to estimate the parameters using the method of moments. So far when the sample size is large at least 20 and with alpha =0.7 and beta=>1.5 then there is no problem. However, If I take n=15 and up to 100 simulations then there is no problem with the code. So I was encouraged to take 1000 simulations and that is when I start having troubles. To answer your second question I am using the first and second moments. I didn't use the variance.

[Rick_SAS]

What is the distribution?

[Rick_SAS · Accepted solution]

Thanks for posting the equations. Instead of solving 

f1(alpha, lambda) = A

f2(alpha, lambda) = B

you should solve

log( f1(alpha, lambda) ) = log(A)

log( f2(alpha, lambda) ) = log(B)

This enables you to use the LGAMMA function instead of the GAMMA function. For more information, see Need to log-transform a distribution? There's a SAS function for that! - The DO Loop

start MOM_Fun(var) global (X);
   X_bar=X[:];
   Xsqr_bar=(X##2)[:];

   alpha = var[1]; lambda = var[2];
   f = j(1, 2, .); /* return a ROW VECTOR */
   f[1] = log(X_bar) - 
          (lgamma(alpha+1) + lgamma(1/lambda+1) - lgamma(alpha+1/lambda+1));
   f[2] = log(Xsqr_bar) - 
           (log(alpha) + lgamma(alpha) + lgamma(2/lambda+1) - lgamma(alpha+2/lambda+1));
   return (f);
finish;

[Salah]

That is brilliant thank you so much