Trapezoidal Rule in PROC IML

Question

Build an IML module tr that calculates an integral (“area under the curve”) using the trapezoidal rule. The module should accept vectors x and y as input parameters and return the value of auc. Show that your module works by applying it as follows:

x = do(-2,5,0.01);         
print "Integral of a over x is" ( tr(x,0#x+1) );         
print "Integral of b over x is" ( tr(x,1#x) );         
print "Integral of c over x is" ( tr(x,x##2) );        
print "Integral of d over x is" ( tr(x,abs(x)) );         
print "Integral of y over x is" ( tr(x,sin(x##2)) );         
print "Integral of z over x is" ( tr(x,log(2+exp(x))) );

Here is my code:

proc iml;
start tr(x, y);
    do i=1 to 2000;
    auc = (sum(x[i]-(x[i-1]))#(y[i]+(y[i-1])))/2;
    return (auc);
    end;
finish;

When I try to run the verification code provided, I get the error message that the subscript is invalid or out of date. What does this mean and how do I fix my code so that it works?

Why are you trying to index a vector `x` of length 701 with values between (and including) 0 and 2000? — whuber, May 09 '14 at 13:51
I am still fiddling with that one, but this definitely works. proc iml; start tr(x,y); n = ncol(x); dx = x[2:n] - x[1:n-1]; ymean = (y[2:n] + y[1:n-1]) / 2; return(dx` * ymean ); finish tr; — user44046, May 09 '14 at 14:05

score 0 · Answer 1 · answered May 09 '14 at 13:45

0

The x only has one row, transpose your x by adding `.

answered May 09 '14 at 13:45

score 0 · Answer 2 · edited May 23 '17 at 10:26

There are a couple of issues here. First off, you don't need to do looping in side your module if the arguments are being passed as vectors. Second, you're indexing i from 1 but you use i-1, which evaluates to 0 when i=1. This causes a subscript error since IML vectors are indexed from 1, not 0 like in C. You'll also get a subscript error whenever x and y have length greater than 2000 since you're referencing elements of the vectors which don't exist.

This is a succinct way of getting what you need:

proc iml;
    start tr(x, y);
        i = 2:ncol(x);
        return( (x[i] - x[i-1])` * (y[i] + y[i-1]) / 2 );
    finish;

    <verification code goes here>
quit;

This is a duplicate of this question, but there are no upvoted or accepted answers to that question as of this posting so I'm unable to mark this as a duplicate. It was even posted on the same day.

Also, this smells like homework. If it is homework, please mention so in your post.

Trapezoidal Rule in PROC IML

2 Answers2