topic ERROR in using NLPNRA and PROC IML for MLE of custom function in SAS/IML Software and Matrix Computations

ERROR in using NLPNRA and PROC IML for MLE of custom function

roystriker77r — Mon, 29 Apr 2024 17:42:04 GMT

Dear community,

I am trying to follow Rick Wicklin's example of using NLPNRA with PROC IML for maximum likelihood estimation. Unfortunately, I am encountering an error with NLPNRA.

I am trying to use this script to find the MLE of a two parameter, two variable custome probability mass function:

/*Define Density Function*/

data multivariate_data;
call streaminit(5);
%let theta1 = 5;
%Let theta2 = 3;
	array Xi[2];
    do i = 1 to 1000;
        do j=1 to 2;
        Xi[j]= rand('Bernoulli',.5);
        end;
  	N=exp(&theta1*(Xi[1]+Xi[2])+&theta2*Xi[1]*Xi[2]); /*numerator of density function*/
  	D=(1+2*exp(&theta1)+exp(2*&theta1+&theta2)); /*denominator of density function*/
  	P=N/D;
        output;
        
    end;
    drop i j;
run;

proc print data= multivariate_data;
var P;
RUN;


proc iml;
use multivariate_data;
read all var {"Xi1","Xi2"} into x;
close multivariate_data;
/* write the log-likelihood function the distribution of interest */
start LogLik(param) global (x);
   theta1 = param[1];
   theta2 = param[2];
   n = nrow(x);
   return ( theta1*(sum(Xi1)+sum(Xi2)) + theta2*sum(Xi1)*sum(Xi2)-n*log(1+2*exp(theta1)+exp(2*theta1+theta2)) );
finish LogLik;



p = {2 4};/* initial guess for solution */
con = {., .}; /* Upper bounds for parameters (if any) */
optn = {1,  /* find max of function       */
        3}; /* print a lot of output      */
call nlpnra(rc,xres,"LogLik",p,optn,con);

However, after running this code I receive an error code following the nlpnra call that says error in argument BLC. After reviewing the documentation I thought that
BLC was an option argument so I am confused as to why I am receiving this error. Anyt help would be much appreciated.
Thanks!

Error in using NLPNRA and PROC IML for MLE of a custom function

roystriker77r — Mon, 29 Apr 2024 19:45:24 GMT

Dear community, I am trying to follow Rick Wicklin's example of using NLPNRA with PROC IML for maximum likelihood estimation. Unfortunately, I am encountering an error with NLPNRA. I am trying to use this script to find the MLE of a two parameter, two variable custome probability mass function:

/*Define Density Function*/

data multivariate_data;

call streaminit(5);

%let theta1 = 5;

%Let theta3 = 3;

array Xi[2];

do i = 1 to 1000;

do j=1 to 2;

Xi[j]= rand('Bernoulli',.5);

end;

N=exp(&theta1*(Xi[1]+Xi[2])+&theta3*Xi[1]*Xi[2]);

D=(1+2*exp(&theta1)+exp(2*&theta1+&theta3));

P=N/D;

output;

end;

drop i j;

run;

proc print data= multivariate_data;

var P;

RUN;

proc iml;

use multivariate_data;

read all var {"Xi1","Xi2"} into x;

close multivariate_data;

/* write the log-likelihood function the distribution of interest */

start LogLik(param) global (x);

theta1 = param[1];

theta3 = param[2];

n = nrow(x);

return ( theta1*(sum(Xi1)+sum(Xi2)) + theta3*sum(Xi1)*sum(Xi2)-n*log(1+2*exp(theta1)+exp(2*theta1+theta3)) );

finish LogLik;

p = {2 4};/* initial guess for solution */

con = {., .}; /* Upper bounds for parameters (if any) */

optn = {1, /* find max of function */

3}; /* print a lot of output */

call nlpnra(rc,xres,"LogLik",p,optn,con);

However, after running this code I receive an error code following the nlpnra call that says error in argument BLC. After reviewing the documentation I thought that BLC was an option argument so I am confused as to why I am receiving this error. Anyt help would be much appreciated. Thanks!

Re: ERROR in using NLPNRA and PROC IML for MLE of custom function

Rick_SAS — Mon, 29 Apr 2024 20:00:31 GMT

Yes, the BLC parameter is optional, and you opted to specify it by using the CON argument. But the CON argument is not the correct format. If you don't want a constrained optimization, then omit the CON argument. Otherwise, specify it correctly.

I don't know what distribution you are trying to fit, but I want to point out that your input data set only has three unique pieces of information:

When (x1,x2)=(0,0), the value P(x1,x2)=2.259E-6
When (x1,x2)=(1,1), the value P(x1,x2)=0.999
When (x1,x2)=(1,0) or (0,1), the value P(x1,x2)=0.00033

I also do not think you can perform an optimization in this function because it is unbounded. For example, look at the cross section for theta1=3 as you let theta2 vary:

/* visualize the LL when theta1=3 */
theta1 = 3;
theta2 = T( do(0,8,0.1) );
LL = j(nrow(theta2), 1);
do i = 1 to nrow(theta2);
   LL[i] = LogLik( theta1 || theta2[i] );
end;
title "LL(3, theta2)";
call series(theta2, LL) grid={x y};

As you can see, the objective function does not have a maximum. I suggest you check the computations that you used to construct the LogLik function.