-1

我尝试运行以下命令以进行数字积分：

nu = 8; 
psi=-0.2;
lambda = 1;
git = @(u) tpdf((0 - lambda * skewtdis_inverse(u, nu, psi)), nu);
g(t,i) = integral(git,1e-10,1-1e-10,'AbsTol',1e-16);

其中 tpdf 是一个 matlab 函数，而 skewtdis:inverse 如下所示：

function inv = skewtdis_inverse(u, nu, lambda)
% PURPOSE: returns the inverse cdf at u of Hansen's (1994) 'skewed t' distribution

c = gamma((nu+1)/2)/(sqrt(pi*(nu-2))*gamma(nu/2));
a = 4*lambda*c*((nu-2)/(nu-1));
b = sqrt(1 + 3*lambda^2 - a^2);

if (u<(1-lambda)/2);
inv = (1-lambda)/b*sqrt((nu-2)./nu)*tinv(u/(1-lambda),nu)-a/b;
elseif (u>=(1-lambda)/2);
inv = (1+lambda)/b*sqrt((nu-2)./nu).*tinv(0.5+1/(1+lambda)*(u-(1-lambda)/2),nu)-a/b;
end

我得到的是：

skewtdis_inverse 错误（第 6 行） c = gamma((nu+1)/2)/(sqrt(pi*(nu-2))*gamma(nu/2));

在调用“F:\Xyz\skewtdis_inverse.m>skewtdis_inverse”期间未分配输出参数“inv”（可能还有其他参数）。

@(u)tpdf((0-lambda*skewtdis_inverse(u,nu,psi)),nu) 中的错误

积分计算/迭代标量值错误（第 314 行）fx = FUN(t)；

积分计算/vadapt 中的错误（第 133 行）[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

积分计算中的错误（第 76 行）[q,errbnd] = vadapt(@AtoBInvTransform,interval);

积分错误（第 89 行） Q = integralCalc(fun,a,b,opstruct);

但是，如果 i 直接调用 thr 句柄中的函数，则没有问题：

tpdf((0 - lambda * skewtdis_inverse(1e-10, nu, psi)), nu)

答案=

1.4092e-11

tpdf((0 - lambda * skewtdis_inverse(1-1e-10, nu, psi)), nu)

答案=

7.0108e-10

非常感谢您的努力！

4

1 回答 1

0

默认情况下，integral期望函数句柄接受向量输入。在您的代码中，if-statement 会产生复杂性，因为仅当满足 ittrue 的所有元素时，u条件才会评估为。因此，如果u是一个具有大于和小于元素的向量(1-lambda)/2，inv则永远不会被分配。

有两种选择：

  1. 将if-statement 放入for-loop 并遍历u.
  2. 使用逻辑索引进行分配。

对于大元素计数，第二个选项更快，并且在我看来，更清洁：

inv = u; % Allocation

IsBelow =  u <  (1-lambda)/2; % Below the threshold
IsAbove =  ~IsBelow         ; % Above the threshold

inv(IsBelow) = (1-lambda)/b*sqrt((nu-2)./nu)*tinv(u(IsBelow)/(1-lambda),nu)-a/b;
inv(IsAbove) = (1+lambda)/b*sqrt((nu-2)./nu)*tinv(0.5+1/(1+lambda)*(u(IsAbove)-(1-lambda)/2),nu)-a/b;
于 2014-10-23T21:43:40.700 回答