matlab - MATLAB奇怪的错误Gamma函数数值积分

Question

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

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)中的错误

积分计算器/iterateScalarValued中的错误(第314行)
fx = FUN(t);

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

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

积分错误(第89行)
Q =积分计算(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

非常感谢您的努力!

Answer 1

默认情况下，integral期望函数句柄接受 vector 输入。
在您的代码中，if -statement创建了一个复杂的函数，因为仅当true的所有元素都满足时，条件的评估结果才为u 。
因此，如果u是一个 vector ，且其元素都大于和小于(1-lambda)/2，则将永远不会分配inv。

有两种选择:

将if语句放入for -loop中，并遍历u的所有元素。

使用逻辑索引进行分配。

对于大量元素，第二种选择更快，我认为更干净:

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;