an integral puzzle

[runningneptune]

Hi all,

I created a macro and used IML in this macro, then I ran this macro twice trying 2 different sets of parameters. I understand that the integral is NOT integrable so the very last column in the matrix should be missing. That said, the 1st try gives me this expected result without any error msgs, but the 2nd does not even work at all, why?

%macro setup(
n1=,
n2=,
r1=,
r2=,
alpha1=,
beta1=,
alpha2=,
beta2=
);

proc iml;

start grand(u) global(a1,a2,a3,a4,a5,a6);
v=(u##(a1-1)#(1-u)##(a4-a1-1))#(1-u#a5)##(-a2)#(1-u#a6)##(-a3);
return(v);
finish;

start Appell(v_a1,v_a2,v_a3,v_a4,v_a5,v_a6) global(a1,a2,a3,a4,a5,a6);
a1=v_a1;
a2=v_a2;
a3=v_a3;
a4=v_a4;
a5=v_a5;
a6=v_a6;
call quad(z,"grand",{0 1}) peak=0;
f=gamma(a4)/gamma(a1)/gamma(a4-a1)#z;
return(f);
finish;

den=j(41,8,.);
do i=0 to 40;
k=i+1;
den[k,1]=i;
den[k,2]=&beta2+&n2-den[k,1];
den[k,3]=&alpha1+&beta1+&n1+&alpha2+&beta2+&n2-2;
den[k,4]=1-(&alpha2+den[k,1]);
den[k,5]=&alpha1+&r1+&beta2+&n2-den[k,1];
den[k,6]=1;
den[k,7]=1;
den[k,8]=Appell(den[k,2],den[k,3],den[k,4],den[k,5],den[k,6],den[k,7]);
end;

print den;

quit;
%mend;

/* 1st try */

%setup(
n1=40,
n2=40,
r1=4,
r2=14,
alpha1=4,
beta1=36,
alpha2=1,
beta2=1
);

/* 2nd try */

%setup(
n1=40,
n2=40,
r1=10,
r2=20,
alpha1=10,
beta1=30,
alpha2=1,
beta2=1
);

[Rick_SAS]

It is the same problem that I told you about when you posted the question to my blog. When u is near 1, the third term in your integrand [(1-u#a5)##(-a2)] will overflow. You can't raise a small number to a very negative power.

To perform the integral, you must be able to evaluate the integrand at every point in the interval (0,1). Otherwise, you will get an error.

In the second call, when k = 30, the constants are as below. When (1-u#a5)=1/1000, then the third term is (1E3)#120 = 1E360, which exceeds the largest floating-point double value.

a1 = 11;
a2 = 120;
a3 = -30;
a4 = 31;
a5 = 1;
a6 = 1;
/* the exponentials */
e1 = a1 - 1;
e2 = a4-a1-1;
e3 = -a2;
e4 = -a3;
print e1 e2 e3 e4;
/* can we evaluate the integrand when u is close to 1? No! */
u = 0.999;
v=(u##(a1-1) # (1-u)##(a4-a1-1)) # (1-u#a5)##(-a2) # (1-u#a6)##(-a3);

[runningneptune]

Hi Rick,

First, I apologize, I didn't mean to double-post; I actually tried different things before I re-posted here. That said, sorry about the bugging on your blog.

I read through your responses, much appreciated! I agree that this integral is not integrable, and I get that. But what strikes me is that why the first try works fine without any error msgs (although the integral did not work, which is expected 🙂 ) ; but the second one did not even work with error msg saying "Invalid argument ..."?

I know that this integration is not valid, but one try tells me error msg; the other one did not but just gave me missing, which is expected.

Thanks a lot!

[Rick_SAS]

And since you have asked several questions like this, here is some general programming advice. The first two will make your programming and debugging easier:

1. Never try to solve an integral on (a,b) until you have successfully evaluated the integrand on (a,b).

2. Never embed a program in a macro until it runs successfully in open code.

3. Reread the article "Guiding numerical integration: The PEAK= option"

 for tips for choosing the PEAK= option. The article says "do not use an endpoint" and "do not use a location where the integrand is zero." You have set PEAK=0, which violates both conditions.

[Rick_SAS]

Because the first set of parameter values do not cause the integrand to overflow (at the points at which the integrand is evaluated). For example, on the 40th iteration of the loop, the parameter values are as below:

/* 1st try */
a1 = 1;
a2 = 120;
a3 = -40;
a4 = 9;
a5 = 1;
a6 = 1;

/* 2nd try */
/*
a1 = 11;
a2 = 120;
a3 = -30;
a4 = 31;
a5 = 1;
a6 = 1;
*/
e1 = a1 - 1;
e2 = a4-a1-1;
e3 = -a2;
e4 = -a3;
print e1 e2 e3 e4;

u = 0.999;
v=(u##(a1-1) # (1-u)##(a4-a1-1)) # (1-u#a5)##(-a2) # (1-u#a6)##(-a3);

Solving an integral is easier when you know what the integral means and how it was constructed. What does this integral represent? Is it related to a probability distribution? (For example, a beta distribution?) Why are the powers in the integrand so large? Where do the parameter values come from?

In similar problems that I have worked on, the parameters were obtained by fitting a model to data. When the model does not fit the data, the parameter estimates can be extreme and the integral might not converge. In these cases, it is not the integral that requires our attention, it is the modeling process that needs revising.