PROC IML梯形规则

【问题标题】：PROC IML trapezoidal rulePROC IML梯形规则
【发布时间】：2014-05-09 05:13:19
【问题描述】：

proc iml;
start tr(x,y); * create function called tr;

N = nrow(x);
dx = x[2:N] - x[1:N-1];
ymean = (y[2:N] + y[1:N-1]) / 2;
return(dx` * ymean );
finish tr;

x = do(-2,5,0.01);
print "Integral of a over x is" ( tr(x,0#x+1) ); *Answer is 7;

我不断收到（执行）无效下标或下标超出范围。我该如何解决这个问题并得到正确的答案？我试过取出 x[1:N-1] 中的 -1;和 y[1:N-1]，但它给了我错误的答案。是因为我需要假设等距间隔吗？如果是这样，我会怎么做。梯形方程：(x-x0)*(y+y0)/2 或 (xi - xi-1) * (yi + yi-1) / 2。

【问题讨论】：

标签： sas calculus sas-iml

【解决方案1】：

下标问题源于x 被视为行向量而不是列向量，因此nrow(x) = 1。这意味着您已定义 n = 1，因此当索引为 n-1 时，它实际上为 0。这会导致下标错误，因为 SAS IML 向量的索引是从 1 而不是 0。要修复它，请使用 n = ncol(x)。

另外，顺便说一句，您可以根据需要缩短模块 tr：

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

    x = do(-2, 5, 0.01);

    print 'The integral of a over x is' (tr(x, 0#x+1)) 'and should be 7.';
quit;

【讨论】：