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
- Permalink
- Email to a Friend
- Report Inappropriate Content

12-10-2013 09:49 AM

I am running a large number of simulations. Below is a snippet of my code. The following code works but is not efficient as a result of the do loop. Can the QUAD function be vectorized for better efficiency?

Regards,

Raphael

/* Incomplete Gamma Function */

start IGF(a,b);

return( b#gamma(a)#CDF('GAMMA', 1, a, 1/b) );

finish IGF;

/* Derivative of the Incomplete Gamma Function w.r.t. aa */

start DIGF(x) global(aa,bb);

return( log(x)#x##(aa - 1/2)#exp(-x#bb/2) );

finish DIGF;

nu = 3;

seed = 123456;

c = j(15000, 1, seed);

w = 2*uniform(c);

evalnum = J(nrow(w),1); /* pre-allocate matrix */

limits = {0 1};

do k = 1 to nrow(w);

aa = nu; bb = w

call quad(result, "DIGF", limits);

evalnum

end;

evalden = IGF(nu + 1/2,w/2);

r = evalnum/evalden;

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

Posted in reply to R_Fraser

12-10-2013 10:53 AM

The answer to your question is no. The integral is a scalar quantity. As you change the parameters, the shape of the function changes and therefore the integration algorithm is different for each parameter. I suspect that the work required to set up each integral is small compared to the work to compute the integral, so I wouldn't expect much gain if QUAD were vectorized (maybe less than 10%?). So I don't think that the DO loop is costing you much time.

A question about your code: Your comment says "Derivative of the incomplete gamma fnc w.r.t. a"

Are you sure that formula is correct? I don't see where the 2s come from. That is, I wonder why (d/da)(IGF) would generate the terms x##(aa-1/2)#exp(-x#bb/2), rather than x##(a-1)#exp(-t##bb). But I admit that I don't understand what you are attempting.

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

Posted in reply to Rick_SAS

12-10-2013 12:15 PM

Thank you. The incomplete gamma function I am using is \int_{0}^{1} u^{(a + 1/2) - 1} \exp(-u*b) du and not the usual \int_{0}^{x} u^{a - 1} \exp(-u) du.

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

Posted in reply to R_Fraser

12-10-2013 03:13 PM

Or rather \int_{0}^{1} u^{(a + 1/2) - 1} \exp(-u*b/2) du