how to write the Constraints when one of the parameters depends on the random variable

[Salah]

Hello ...

I am trying to solve a system of nonlinear equations. My problem is with one of the constrains that depends on X. My X is a vector with dimensions n x 1. I looked on the website thinking that I may find something I couldn't. I would really appreciate any guidance.

Thank you

[PeterClemmensen]

Did you read this?

https://blogs.sas.com/content/iml/2018/02/28/solve-system-nonlinear-equations-sas.html

[Salah]

Yes I did and from it I learned how to construct the equations. But in the link you provided the constrains are numbers no variables like in my case

[Rick_SAS]

One way to do this is to check the condition in the objective function and add a penalty if the condition is violated.

For example, if you require that X > mu, then compute 

b = ^(X > mu);   /* b=1 if condition true */

penalty = b * 1e6 * (X - mu);

Then add the penalty to the objective function.

[Salah]

Hello Rick,

Thank you for the reply.  Do you mean like this

b = ^(X > mu); 

penalty = b * 1e6 * (X - mu);

con = {1e-6   1e-6   -inf,   .    .   penalty };

But what if X is a vector that doesn't have one value would that still work?

Thank you

[Rick_SAS]

No. THe BLC matrix is only for linear and boundary  constraints on the parameters.

Check the condition in the objective function, which you call 'FUN', and return a large value if the current value of mu does not satisfy the condition on X.

Now that I look at your implementation, I see that you can't even evaluate the objective function if X < mu, because the objective function has terms like log(X - mu).. What are you trying to do? It sort of looks like a maximum likelihood optimization (sine you simulate X as a sample from the Weibull distribution), but an MLE compuation would use a scalar objective function.

Here is a link to some resources:
- MLE estimation example in SAS/IML

- How to use the LOGPDF for common distributions 

- 10 tips to follow before you run an optimization

[Salah]

Thank you Rick for your help and for your time.

My problem is not to find the MLEs. I have 3 equations that I need to solve which I defined in f(1)-f(3). I have to solve other system of equations that are almost have same nature ( X-Mu)

I am afraid that I didn't quit understand what you said....I am copying that part of my code below. Kindly ask you to look at it and let me know if what I did is correct and acceptable.  Thank you.

start fct(var) global (X,n);

Alpha = var[1]; Mu = var[2]; Theta = var[3];

f = j(1, 3, .); 

sum1=0;sum2=0;sum3=0;

b=J(n,1,.);
b=X-Mu;

do i=1 to n;

if b[i]<=0 then b[i]=0.00001;
if b[i]>0 then b[i]=X[i]-mu;
end;

do j=1 to n;

sum1=sum1+log(b[j]);
sum2=sum2+ ( b[j])**Alpha;
sum3=sum3+ log(b[j]) *( b[j] )**Alpha;

end;

f[1] = ( log(theta)+digamma(1.0))/alpha - sum1/n;

f[2] = sum2/n - theta;

f[3] = ( log(theta)+digamma(1.0)+1)/alpha + sum3/(n*theta) ;

return (f);
finish;

[Rick_SAS]

Well, it's hard for me to determine if your program is correct or not because I don't understand what statistical problem you are trying to solve. However, my suggestion was for you to use

b=X-Mu;
if any(b<0) then 
   return( {1e20, 1e20, 1e20} );

in the FCT function immediately after you define the b vector.

Even if you are not doing MLE, I suggest you study the IML programs in the links I sent. There are several suggestions and best practices in the examples that I think would improve your program.

[Salah]

I am sorry if I wasn't clear about my problem. I am trying to use the principle of maximum entropy to estimate the parameters of the Weibull distribution.

To do this problem, I have to first  to construction the Lagrange multiplier (lambda0, lambda1 and lambda2). Then I derive the entropy function and find a connections/relation between the parameters and these lambdas.  After that I find the derivatives of  the entropy function with respect to lambda0, lambda1 and lambda2. Finally I try to use the data and the relationship b/w these lambdas  and the parameters which I get in terms of  3 nonlinear equations that I need to solve simultaneously.

What I did in my code is to simulate data from Weibull then wrote these equations in subroutine and try to solve them. I used your recommendations unfortunately I keep getting the following message

ERROR: HYQUAN Optimization cannot be completed.
ERROR: HYQUAN needs more than 200 iterations or 500 function calls.

I tried to use tc={10000 15000}; and define it inside the call function (see attached code), this time I received this warning message

ERROR: HYQUAN Optimization cannot be completed.
WARNING: Optimization routine cannot improve the function value.

Thank you,


 

[Rick_SAS]

The following program is based on your code. I made a few efficiency improvements to eliminate unnecessary loops.

*** Solve a system of nonlinear equations in SAS/IML  **;
*** reference: https://blogs.sas.com/content/iml/2018/02/28/solve-system-nonlinear-equations-sas.html;

proc iml;
alpha=2; mu=4;  theta=1;  /*Assumed real values */
n=20; seed=2; k=100;

start fct(var) global (X,n);
   Alpha = var[1]; Mu = var[2]; Theta = var[3];
   f = j(1, 3, .);         /* return a ROW VECTOR */
   sum1=0;sum2=0;sum3=0;
    
   b=X-Mu;
   if any(b<0) then 
      return( {1e20, 1e20, 1e20} );

   sum1 = sum( log(b) );
   sum2 = sum( b##Alpha );
   sum3 = sum( log(b) # b##Alpha );

   f[1] = ( log(theta)+digamma(1.0))/alpha - sum1/n;
   f[2] = sum2/n - theta;
   f[3] = ( log(theta)+digamma(1.0)+1)/alpha + sum3/(n*theta) ;
   return (f);
finish;
 
Alpha1=J(k,1,.);mu1=J(k,1,.);theta1=J(k,1,.);

/* These are the parameters of the Weibull(theta,alpha,mu):   Theta>0      Alpha>0      Mu<X  */   
con = {0.01    0.01  0.01,  
       200     200   200}; 
     
x0 = {1.5 3 2};           /* initial guess */
optn = {3, 0};               
tc={10000 15000};

Do N1=1 to k;
   U = Uniform( J(n,1,seed) );
   X = ( - Theta* Log(1 - U ))##(1/Alpha) + Mu; 

   call nlphqn(rc, Soln, "fct", x0, optn, con,tc) ;
   if rc > 0 then 
      Alpha1[N1] = Soln[1]; Mu1[N1] = Soln[2]; Theta1[N1] = Soln[3];
   *print Soln;
End;

alpha_hat=Alpha1[:];
mu_hat=Mu1[:] ;
theta_hat=Theta1[:];

*print Alpha1 Mu1 Theta1;
print alpha_hat  mu_hat theta_hat ;

[Salah]

Thank you so much Rick.

I highly appreciate that  you were so patient with me, God bless you.

Do I have to worry about the message below??

ERROR: HYQUAN Optimization cannot be completed.
WARNING: Optimization routine cannot improve the function value.

Does the above error message has anything to do with initial values?  I looked at the link that you provided one of them is using SAS/studio to see the 3D graph which I don't have. Is there any other alternatives to picture the pdf and the MLEs?

Thank you again for everything you do for us.

[Rick_SAS]

Are you seeing that error when you run EXACTLY the code I posted? Because I do not get that error when I run SAS/IML 15.1.

Regardless, if you see that message a few times among the 100 simulated samples, it means that the optimization did not converge for certain samples. This can happen if your initial guess is not very good or if you are using very small simulated samples. (I classify N=20 as small for a three-parameter family.) For small simulated samples, you might get "unlucky" and the sample does not look very much like a sample from the distribution with the specified parameters...

There are two ways to handle the situation. 

1) Make a better initial guess 

2) Report the results of the simulation by saying, "Of the k simulated samples, three did not converge. The results for the remaining samples are shown in Table 3...."

You shouldn't get that error (or not very often) for larger samples like N=100 or more.

[Salah]

Yes i got the error message running the exact same code ( I am using SAS 9.4). I did try to work on the initial values by using one of the codes that you posted on this link. I succeeded for some values and I didn't for others.

I found this message 

ERROR: Expecting page 3, got page -1 instead.
ERROR: Page validation error while reading WORK.'SASTMP-000000001'n.ITEMSTOR.
ERROR: File WORK.'SASTMP-000000001'n.ITEMSTOR is damaged. I/O processing did not complete.

At the beginning I thought that this message is due to getting negative likelihood function but even with +ve log-likelihood I get this message. 

I am not sure if I should omit the values that cause these error messages but in this case I won't cover a large spectrum of values?

[Rick_SAS]

PROC IML does not write itemstores, so that error message is coming from somewhere else. How are you running SAS? Are you using SAS Studio? Enterprise Guide? Are you using the SAS University Edition (UE free version)? The message says that the item store is corrupted, so maybe there was an installation issue.  I don't know.

I did not get any errors on SAS 9.4 on Windows. Are you running Linux (or on the SAS UE virtual machine?)

Again, for small simulated samples, you might occasionally get these nonconvergence issues. It will be worse in some regions of parameter space than in others. That is part of what a simulation study tells you, and you should report the results accordingly. FOr example, "for samples of size 20, the optimization only converges 80% of the time when alpha is smaller than [or larger than] ...."

[Rick_SAS]

Is this problem resolved or are there still problems? If it is resolved, please mark it as Solved.

I wrote down some thoughts and tips about how to define constraints when the parameters in an optimization are constrained by the data and have a restricted domain. See the article "Two tips for optimizing a function that has a restricted domain."