How to perform integral in SAS that can not be evaluated in closed form

Accepted Solution Solved
Reply
Super Contributor
Posts: 297
Accepted Solution

How to perform integral in SAS that can not be evaluated in closed form

How do I perform integration in SAS that can not be evaluated in closed form


Accepted Solutions
Solution
‎12-19-2013 09:46 AM
Respected Advisor
Posts: 2,655

Re: How to perform integral in SAS that can not be evaluated in closed form

If you have access to SAS/ETS, look at PROC MODEL, and in particular at Example 19.7 Spring and Damper Continuous System and Example 19.10 Systems of Differential Equations.  These may provide a framework for your particular question.

Steve Denham

View solution in original post


All Replies
Solution
‎12-19-2013 09:46 AM
Respected Advisor
Posts: 2,655

Re: How to perform integral in SAS that can not be evaluated in closed form

If you have access to SAS/ETS, look at PROC MODEL, and in particular at Example 19.7 Spring and Damper Continuous System and Example 19.10 Systems of Differential Equations.  These may provide a framework for your particular question.

Steve Denham

Super Contributor
Posts: 297

Re: How to perform integral in SAS that can not be evaluated in closed form

Hi my proc model has no example 19.7, it is just chapter 14

Respected Advisor
Posts: 2,655

Re: How to perform integral in SAS that can not be evaluated in closed form

The subsections will have the same titles.  This capability has been in PROC MODEL since 9.1.3, I believe.

Steve Denham

Respected Advisor
Posts: 4,646

Re: How to perform integral in SAS that can not be evaluated in closed form

Look here and here. - PG

PG
SAS Super FREQ
Posts: 3,477

Re: How to perform integral in SAS that can not be evaluated in closed form

You don't give the dimension of the domain. For one- or two-dimensional integration, use the QUAD function in SAS/IML: How to numerically integrate a function in SAS - The DO Loop

New Contributor
Posts: 3

Re: How to perform integral in SAS that can not be evaluated in closed form

A lot of folks don't have ETS or IML. Here is a routine for numerical integration in base SAS. It could stand some improvements and efficiencies - as many of you know, clever is not my long suit - - but it's understandable, teachable and in base SAS. I wrote it to integrate over spectrum lines in astrophysics which have no closed or even parametric form. The non-parametric f(x) from a to b = [ f(a) + f(a+Δ) + f(a+2Δ) + f(a+3Δ) + ... + f(b) ] * Δ. That is,

1. Break up the interval over which you wish to integrate into tiny intervals. Remember the width of the tiny intervals measured along the x-axis - that's Δ. Often, the data is already in this form: a long list of x-y pairs. If not but you do have a the formula f(x), create just such a set of ordered pairs:

%macro partition(a,b,n);

data partition;

   do i = 0 to &n.;

      x = &a. + i * ((&b. - &a.) / &n.);

      y = f(x);

      output;

   end;

   keep x y;
run;
%mend partition;


In this macro, a is the lower limit of integration, b is the upper limit and n is the number of tiny divisions for the numeric integration. I often use 500; use what gives the level of accuracy you need. This is, of course, the simple rectangle rule written into base SAS. choosing a large enough n should give the accuracy you need.

f(x) is the function on x. NB: this code will not run as written - the f(x) you want must be typed in first. This will produce a set of order pairs x and y=f(x) over the interval from x=a to x=b, if you do not already have this at the start.

2. Once you have such a set of ordered pairs, add together the y values using a PROC SUMMARY and multiply the result by the width of the interval you chose in part 1 (or is native to the data):

PROC SUMMARY data=partition;
   var y;

   output sum=y_sum out=integration;

run;

data integration;

   set integration;
   delta = (b - a) / n;

   integration = y_sum * delta;

run;

where, once again, a and b are the limits of integration and n is the numbmer of intervals. you might want to put all this inside a single macro:

%macro partition(a,b,n);

data partition;

   do i = 0 to &n.;

      x = &a. + i * ((&b. - &a.) / &n.);

      y = f(x);

      output;

   end;

   keep x y;
run;

PROC SUMMARY data=partition;
   var y;

   output sum=y_sum out=integration;

run;

data integration;

   set integration;
   delta = (&b. - &a.) / &n.;

   integration = y_sum * delta;

run;

%mend partition;

Respected Advisor
Posts: 4,646

Re: How to perform integral in SAS that can not be evaluated in closed form

Well explained solution! The proposed integration formula will however give slightly biased results (you are adding n+1 terms and dividing by n). The call to proc summary, simply to calculate a sum, is not necessary either. Everything can be done in one step and without bias with the trapezoid rule. Something like:

%macro integrate(f,a,b,n);
data _null_;
sum = &f(&a);
do i = 1 to &n - 1;
     sum + 2*&f(((&n-i)*&a + i*&b)/&n);
     end;
sum + &f(&b);
integral = sum * abs(&b - &a) / (2*&n);
put "Integral of &f, from &a to &b, using &n steps = " integral;
run;
%mend;

%integrate(sin,0,3.14159265358979,100);

%integrate(sin,0,3.14159265358979,10000);

PG

PG
Super Contributor
Posts: 297

Re: How to perform integral in SAS that can not be evaluated in closed form

Thanks David

New Contributor
Posts: 3

Re: How to perform integral in SAS that can not be evaluated in closed form

PGStats - good catch! It should be n-1 - my mistake when I typed it into the web site. And, of course, use of the trapezoidal rule in your post is better.

Super Contributor
Posts: 297

Re: How to perform integral in SAS that can not be evaluated in closed form

THanks PGStats, thanks a lot

Super Contributor
Posts: 297

Re: How to perform integral in SAS that can not be evaluated in closed form

I guess on of you  could help me with this. Thanks for all the help:

Consider the following .:

Consider the following .

C= know value

U is Uniform (0,1)

F(x)=1-S(x) is known

a=known value

k=known value

ht-= is known.

The only unknown is X.

I wish to write a SAS code that find X such that the right hand sight is equal the left hand side numerically.

I can not make X the subject but can find a numerical solution.

How can I find the Values of X s. that

C*U- [F(x) *exp(-(k-a*X)**2 - (k - a**2 * X) exp(-(k -a*X)**2/2t)*F(x)- ht]=0

☑ This topic is SOLVED.

Need further help from the community? Please ask a new question.

Discussion stats
  • 11 replies
  • 568 views
  • 2 likes
  • 5 in conversation