topic QUAD Function in SAS/IML Software and Matrix Computations

QUAD Function

R_Fraser — Tue, 10 Dec 2013 14:49:00 GMT

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 = result;

end;

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

r = evalnum/evalden;

Re: QUAD Function

Rick_SAS — Tue, 10 Dec 2013 15:53:52 GMT

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.

Re: QUAD Function

R_Fraser — Tue, 10 Dec 2013 17:15:03 GMT

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.

Re: QUAD Function

R_Fraser — Tue, 10 Dec 2013 20:13:34 GMT

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