Statistical programming, matrix languages, and more

Product of integrals and maximization

Reply
Contributor
Posts: 41

Product of integrals and maximization


Hi all,

I have posted parts of the problem before, but am reposting here to discuss the entire problem.

I want to maximize the product of 130 integrals, where the integrans is a complicated binomial density, to get an MLE estimate for a parameter (rho).

Please give me advice on the code structure below. I have split it into 3 steps: define the integrans, define a function to evaluate the integral(using QUAD in IML) and take the product, and finally call an NLP function.

Thanks for any help.

Thanks so much for your reply.

I am posting the entire code for my problem: Obtain the MLE of rho, by maximizing a product of 130 integrals. The integrand is a binomial density function, and one argument to this density function is rho.

There are 3 steps, first to define a function to evaluate the integrand, second to evaluate the integral and the product of all integrals, and finally the maximization.

Below rho , N,D,val,ni,di are all Global variables.

I would like to know if the logic below is fine, or is there anything else I need to consider.

proc iml;


/* Step 1. */

/*** Data ***/

/*** N, D are the input data ***/
/*** There are 130 observations of N and D ***/
/*** There are 130 intergals to be evaluated***/
/*** (N,D) will be the input to the i-th integral ***/

N = j(130,1,1);             /*** N is a constant array of 130 values of n_t.Initialized to 1's here  ***/
D = j(130,1,1);             /*** D is a constant array of 130 values of d_t.Initialized to 1's here ***/
ni=0;
di=0;

result=0;
gamma = 0.3 ;          /*** gamma is a constant. ***/
val = j(130,1,1);         /***column vector initialized to 1's to hold reslt of the integrations using QUAD ***/
iter = 0;                    /*** Iteration counter initialized to 0   ***/

/* Define the integrand, which will be called by QUAD subroutine.*/
/*** The inputs ni and di hold the i-th element of N and D ***/
/*** The function Func will be called  by Prod_Of_integral 130 times to evaluate the respective integrals***/


start Func(y)  GLOBAL(gamma,rho,ni, di,N,D);
      /***Note: if iter = i then ni = N; di=D***/
              p = CDF("Normal",(gamma - {sqrt(rho)*quantile("Normal", y)} ) / sqrt(1-rho) );
              /***Handle p if it is outside [0,1]...***/
              if p >= 1 then p = 1;
              if p <= 0  then p =0.01;
              f = pdf("Binomial", di, p, ni);
    return( f );

finish;


/* Step 2. */
/* Define a function that computes the product of the integrals.*/

start Prod_Of_Integral(parms) GLOBAL(rho, val,iter,ni,di,N,D);
/*the optimization parameter is:   rho */
        rho = parms[1];
        yinterval = 0 || 1;  /* evaluate integral on [0, 1] */
        do i = 1 to 130;
             /*** reset  iter    ***/
             /*** reset ni and di ***/
            iter = i; ni = N; di = D;

            call quad(result, "Func", yinterval); /* evaluate integral on [0, 1] */

             /*** update result***/
             val =result;
      end;

         return(prod(val));
finish;

/* Step 3. */


/*Maximize the product of the integrals to obtain MLE for rho using nlpnra */

/* Set up linear constraint 0 < rho < 1 */
con = {0 ,1}

/* initial guess for solution for rho */
rho_0 = {0.5};
/* find maximum of function */
opt = {1,2} ;

call nlpnra(rc, result, "Prod_Of_Integral", rho_0, opt) blc=con ;

Thanks again for your advice.

Ask a Question
Discussion stats
  • 0 replies
  • 366 views
  • 0 likes
  • 1 in conversation