turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

Find a Community

- Home
- /
- SAS Programming
- /
- General Programming
- /
- How to perform integral in SAS that can not be eva...

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

12-19-2013 03:32 AM

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

Accepted Solutions

Solution

12-19-2013
09:46 AM

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

12-19-2013 09:46 AM

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

All Replies

Solution

12-19-2013
09:46 AM

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

12-19-2013 09:46 AM

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

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

12-20-2013 01:24 PM

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

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

12-20-2013 01:40 PM

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

Steve Denham

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

12-20-2013 01:43 PM

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

12-26-2013 07:47 AM

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

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

12-27-2013 03:55 PM

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;

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

12-27-2013 11:08 PM

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

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

12-30-2013 12:38 AM

Thanks David

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

12-28-2013 05:44 PM

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.

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

12-30-2013 12:39 AM

THanks PGStats, thanks a lot

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

12-30-2013 12:51 AM

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**