sas - SAS：使用 IML 处理矩阵

Question

所以我试图计算这个公式，但结果很奇怪。元素非常大，所以我不确定我哪里出错了。我附上了公式的照片：

这是我的代码：

*calculating mu_sum and sigma_sum;
T_hat=180;
mu_sum_first_part={0,0,0,0};
mu_sum_second_part={0,0,0,0};
mu_sum={0,0,0,0};

*calculating mu_sum;
do i = 0 to T_hat;
    term=(T_hat - i)*(B0**i)*a;
    mu_sum_first_part = mu_sum_first_part + term;
end;
do i=1 to T_hat; 
    term =B0**i;
    mu_sum_second_part = mu_sum_second_part + term;
end;
mu_sum = mu_sum_first_part + mu_sum_second_part*zt;
print mu_sum;

*calculating sigma_sum;
term=I(4);
sigma_sum=sigma;
do j=1 to T_hat;
    term = term + B0**j;
    sigma_sum = sigma_sum + (term*sigma*(term`));
end;
print sigma_sum;

我知道这很长，但请帮忙！！

Answer 1

跳出来的第一件事是你的循环第一个术语mu有 1 太多：

do i = 0 to T_hat;
    term=(T_hat - i)*(B0**i)*a;
    mu_sum_first_part = mu_sum_first_part + term;
end;

应该：

do i = 0 to T_hat-1;
    term=(T_hat - i)*(B0**i)*a;
    mu_sum_first_part = mu_sum_first_part + term;
end;

Answer 2

您的程序在数学上没有任何问题。当您将矩阵提高到 180 次方时，您应该不会对看到非常大或非常小的值感到惊讶。例如，如果让 B0 = { 0 1 0 0, 0 0 1 0, 0 0 0 1, 0 1 1 1 }; 那么 B0**T 的元素是 O( 1E47 )。如果将 B0 除以 2 并将结果提高到 180 次方，则元素为 O(1E-8)。

据推测，这些公式适用于具有特殊结构的矩阵 B0，例如 ||B0**n|| --> 0 作为 n --> 无穷大。否则幂级数不会收敛。我建议您仔细检查您使用的 B0 是否满足参考假设。

您没有询问效率，但明智的做法是使用SAS/IML 中的霍纳方法计算截断幂级数，而不是明确形成 B0 的幂。