BookmarkSubscribeRSS Feed
runningneptune
Fluorite | Level 6

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
);

4 REPLIES 4
Rick_SAS
SAS Super FREQ

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
Fluorite | Level 6

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
SAS Super FREQ

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
SAS Super FREQ

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.

sas-innovate-2024.png

 

Secure your spot at the must-attend AI and analytics event of 2024: SAS Innovate 2024! Get ready for a jam-packed agenda featuring workshops, super demos, breakout sessions, roundtables, inspiring keynotes and incredible networking events.

 

Register by March 1 to snag the Early Bird rate of just $695! Don't miss out on this exclusive offer. 

 

Register now!

Multiple Linear Regression in SAS

Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.

Find more tutorials on the SAS Users YouTube channel.

From The DO Loop
Want more? Visit our blog for more articles like these.
Discussion stats
  • 4 replies
  • 5866 views
  • 0 likes
  • 2 in conversation