QUAD Function

[R_Fraser]

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;

[Rick_SAS]

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.

[R_Fraser]

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.

[R_Fraser]

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