topic Re: Product of functions/Variable scope in SAS/IML Software and Matrix Computations

Product of functions/Variable scope

nstdt — Tue, 15 Jul 2014 13:27:20 GMT

Hi all,

I need to get the product of 100 integral functions (computed using the quad funciton). Here is some high-level code below.

/*** Variables that will be accessed by other functions including QUAD to store the product of itnergrals ***/

/*** val is an array, each of whose elements will hold the value of one integral ***/

val = j(100,1,1);

/*** ni,di, N,D will be global parameters to the function defining the integrand ***/

ni=0;

di=0;

N=j(100,1,1);

D=j(100,1,1);

gamma=k; /*** This is a constant ***/

/* Define the integrand, which will be called by QUAD subroutine.*/
start Func(y) GLOBAL(gamma,rho,ni, di,n,d);
f = ...;

return( f );

finish;

/* Calculate the Product of integrals using Quad */

start Prod_Of_Integral(parms) GLOBAL(rho, val,count,ni,di,N,D);

/*** Call QUAD a 100 tiems to evaluate each of 100 integrals ***/

/*** Results should be i val ***/

do i = 1 to 100;

// Do stuff

// Reset ni and di to the ith value of N and D

ni = N; di = D;

// QUAD below stores its resultant integral value in val ***/

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

end;

do j = 1 to 100;

prod = prod * val;

end;

return( prod );

finish;

My question is: will val be accessible inside QUAD, and will prod hold the resultant product of all the integrals?

Is it possible to reset ni and di to the i-th value in N and D (during the i-th iteration of the DO loop inside Prod_Of_Integral?

I need val to be a global variable, and similarly the parameters ni,di, N,D.

How do I declare the scope of my variables/ structure the functions to calcualte the product of the integrals?

Thanks for any advice and help.

Re: Product of functions/Variable scope

Hutch_sas — Tue, 15 Jul 2014 14:04:06 GMT

You need to modify the line:

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

val will not be replaced by the quad call. do

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

val = result;

Re: Product of functions/Variable scope

Rick_SAS — Tue, 15 Jul 2014 14:18:43 GMT

...and then delete the last DO loop and just say

return( prod(val) );

However, if the integral values are large, you'll be facing an overflow. For example, if each integral has the value 2, then the product will be 2##100. It is numerically safer to work on the log scale and return log(prod(val) = sum(log(val)).

The expression sum(log(val)) won't overflow.

Nothing in what you've said requires val, count, N, or D need to be global variables, but maybe there are additional reasons that are unstated.

Re: Product of functions/Variable scope

nstdt — Tue, 15 Jul 2014 17:09:18 GMT

Thanks for your reply again.

You said that val, count, N, D don't need to be global variables in this problem(I have presented all the necessary information; though this problem is really a continuation of my earlier post about Maximizing Integrals to which you also replied. The third piece here is the Maximization done by calling NLPNRA, where the variable rho will be a global parameter).

Inside the Prod_Of_Integral function, I reset the value of ni and di to equal N and D respectively. If they are not global, how could I access their values from inside a function?

N and D are may 100-length data vectors in the "main" scope, but inside the function Func, I need the only the i-th elements of N and D to be processed in the i-th iteration.(From my earlier post, f is a binomial density. and ni and di are the parameters. I have 100 observations of ni and di).

The idea is that, at the i-th iteration of the DO loop in Prod_Of_Integral, the corresponding i-th call to QUAD will evaluate an integral with integrand f(ni,di) = Bin(ni,di), where f(ni,di) is computed by Func. I need to pass all the pairs (ni,di) to Func successively, before each result of Func is integrated by QUAD.

Then, the result of the call to QUAD(the result of the integration) will be stored in variable result(as per Hutch above) and then stored in val.

Then I take the product over all elements of val, giving the product of the integrals, which is then maximized, when I call NLPNRA:

rho_0={0.5};

call NLPNRA(rc,res,"Prod_Of_Integral", rho_0);

So,in order to access the values of ni and di outside MAIN scope, as above,don't I need them to be global?

Thanks again for your advice and insights!

Re: Product of functions/Variable scope

Rick_SAS — Tue, 15 Jul 2014 17:17:24 GMT

Yes, N and D should be global if Prod_Of_Integral will be an objective function. Sorry for the confusion.

Re: Product of functions/Variable scope

nstdt — Wed, 16 Jul 2014 14:03:40 GMT

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 integrans, 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. */

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

Re: Product of functions/Variable scope

nstdt — Wed, 16 Jul 2014 14:04:23 GMT

Thanks so much for your reply.

There are 3 steps, first to integrans, 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. */

/*** 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;

finish;

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

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.