turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

Find a Community

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

07-13-2014 01:46 PM

Hi all,

The following likelihood function L which is an integral needs to be maximized. I think the integral has to be numerically computed.

L = int { k1 * p^d * (1-p)^ (n-d) dy } , from y=0 to y=1, where p is a function of y, and k1,n,d are constants, k1= n choose d

p(y) = ...

I need to maximize L to obtain MLE estimates for gamma and rho. (Or another option is to fix gamma as a constant too, and estimate rho). Thus gamma and rho are the unknown quantities of interest.

I am having trouble setting up the integral to be computed, and then maximizing it.

Specifically, I am not sure of the syntax to specify the dependence of the integrand on gamma and rho, in a way that these are identified as the parameters to be optimized.

I read Rick Wicklin's blog **Optimizing a function of an integral, **

where he uses NLPNRA to maximize an integral. Can NLPNRA be applied here too? Or does the integral L require a different treatment?

Here is some IML pseudocode.

proc iml;

/* Define the integrand, which will be called by QUAD subroutine.*/

n = n1 /** n1 is a constant**/

d = d1 /** d1 is a const **/

start Func(gamma,rho) GLOBAL n d;

p = ...;

f = ...;

return( f );

finish;

/* Define a function that computes the integral.*/

start Integral(parms);

/* the parameters are: gamma = parms[1]; rho = parms[2]; */

yinterval = 0 || 1;

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

return( val );

finish;

/* Set up linear constraint 0 < rho < 1 */

con = {. 0 , /* no constraint on lower bound of gamma, 0 < rho */

. 1 } /* no constraint on upper bound of gamma, rho < 1 */

/*Maximimize using nlpnra */

p = {-1 1}; /* initial guess for solution */

opt = {1, /* find maximum of function */

2}; /* print a little bit of output */

call nlpnra(rc, result, "Integral", p, opt) blc=con j;

Thank you for any help and advice.

Accepted Solutions

Solution

07-14-2014
11:29 AM

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

07-14-2014 11:29 AM

All Replies

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

07-14-2014 09:51 AM

to Here are some thoughts:

1) FUNC needs to be a function of y. The gamma and rho parameters go into the GLOBAL clause.

2) NORMCDF(...)-->CDF("Normal",...) and INVNORMCDF(...) --> quantile("Normal,...)

3) For some parameter values, p (in FUNC) will be outside the interval [0,1] because you are dividing the CDF by sqrt(1-rho).

Then you will get an error when you evaluate f. You have to modify FUNC to handle the case where p is not in [0,1].

[WRONG..IGNORE THIS

And the most important observation is this:

4) The integrand f is the PDF of the binomial distribution. It is equivalent to PDF("Binomial",d,p,n); Consequently, you do not need to call QUAD to evaluate this integral. The integral of a PDF is the CDF, so you can use CDF("Binomial",...), to evaluate the integral, as shown in this blog post: Optimizing a function that evaluates an integral - The DO Loop

]

Message was edited by: Rick Wicklin IGNORE the last point, which is wrong. I just realilzed that p is a function of y, not a constant.

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

07-14-2014 10:29 AM

Thanks for the reply. Your blog is really helpful.

On (1), I am struggling to express FUNC with y as parameter. What do I do about n and d? How are those passes into FUNC?

Will they also go into the GLOBAL caluse? I somehow need the NLPNRA to recognize that it is gamma and rho that need to be optimized.

Is this right:

start Func(y) GLOBAL gamma,rho;

p = ...;

/*** Do Something to handle case when p is outside [0,1].... ***/

f =...;

return( f );

finish;

Thanks again for any help and advice.

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

07-14-2014 10:35 AM

start Func(y) GLOBAL(gamma,rho,n, d);

p = ...;

/***Do Something to handle case when p is outside [0,1].... ***/

f = pdf("Binomial", d, p, n);

return( f );

finish;

Solution

07-14-2014
11:29 AM

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

07-14-2014 11:29 AM

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

07-14-2014 03:29 PM

Hi again,

Thanks for your replies. I still have some questions, more on the implementation though.

First, it turns out gamma is a constant, and I only need to optimize for the parameter rho.

In this case, it seems like I can't avoid calling quad. And my confusion is that when I write the Integral function following your blog, then I need to specify rho as the parameter that changes. Howevr, the function Func needs the variable of integration "y" as an input, where y is obtained from the INVNORMCFD and goes from 0 to 1.

How can I simultaneously retain y as the variable in Func but rho as the parameter in Integral ?

Here is my revised code below. My confusion is in Step 2. Should I make it Integral(parm), that is, make the Integral function take an input argument parm which will map to rho, or should rho be outside with the GLOBAL statement a below? Does Intergal need an input argument here at all?

I am confused about how to differentiate between y and rho.

I think that NLPNRA will call Integral with parameter rho, and Integral executes QUAD. Then QUAD calls FUNC. But how do I tell QUAD that though y is the input to FUNC, rho is the parameter? In defining Integral, should I use: start Integral GLOBAL(rho);

or should it be: start Integral(y) GLOBAL(rho) ?

proc iml;

/* Step 1. */

/* Define the integrand, which will be called by QUAD subroutine.*/

n = N /*** N is a constant. ***/

d = D /*** d is a constant. ***/

gamma = k /*** gamma is a constant. ***/

start Func(y) GLOBAL(gamma,rho,n, d);

p = CDF("Normal",(gamma - {sqrt(rho)*quantile("Normal", y)} ) / sqrt(1-rho) );

/***Cap and Floor to handle case when p is outside [0,1].... ***/

If p > 1 then p =1;

If p < = then p =0.01;

f = pdf("Binomial", d, p, n);

return( f );

finish;

/* Step 2. */

/* Define a function that computes the integral.*/

start Integral GLOBAL(rho); /**** should this be Integral(y) GLOBAL(rho) ????????? ***/

/* the optimization parameter is: rho */

yinterval = 0 || 1;

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

return( val );

finish;

/* Set up linear constraint 0 < rho < 1 */

con = {0 ,1}

/*Maximimize using nlpnra */

rho_0 = {0.5}; /* initial guess for solution */

opt = {1,2} /* find maximum of function */

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

Thanks again very much.

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

07-14-2014 03:56 PM

To add a further layer of complexity, in the full problem, n and d are both arrays of length say, 100. And I really need to maximize the Sum of log of the Integrals of Func(the sum is over each of the 100 n and d pairs, rather than just one n1 and d1 as I have used here).

If you have any further hints, please do let me know, I really appreciate it.

Thanks again very much!

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

07-14-2014 04:11 PM

It is a little confusing because there are two functions: one for the integral and one for the NLP function. Try this:

start Func(y) GLOBAL(gamma,rho,n, d);

p = ...;

f = pdf("Binomial", d, p, n);

return( f );

finish;

/* Define a function that computes the integral.*/

start Integral(parm) GLOBAL(rho);

rho = parm[1];

yinterval = 0 || 1;

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

return( val );

finish;

In 'Func', y is the "dummy variable" that is used for integration. The QUAD function will evaluate 'Func' at many points in order to compute the integral of Func(y) dy on [0,1].

In 'Integral', the NLP routine will optimize the parm variable. At each iteration, the global variable rho is updated with the value of parm before the QUAD routine is called. This syntax should help you keep things straight.