Search code examples
matlaberror-handlinggamma-function

MATLAB strange error Gamma function numerical integration

i try to run the following in order to integrate numerically:

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);

where tpdf is a matlab function and skewtdis:inverse looks like this:

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

What i get out is:

Error in skewtdis_inverse (line 6) c = gamma((nu+1)/2)/(sqrt(pi*(nu-2))*gamma(nu/2));

Output argument "inv" (and maybe others) not assigned during call to "F:\Xyz\skewtdis_inverse.m>skewtdis_inverse".

Error in @(u)tpdf((0-lambda*skewtdis_inverse(u,nu,psi)),nu)

Error in integralCalc/iterateScalarValued (line 314) fx = FUN(t);

Error in integralCalc/vadapt (line 133) [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 76) [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 89) Q = integralCalc(fun,a,b,opstruct);

If i , however call the function in thr handle directly there are no Problems:

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

ans =

1.4092e-11

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

ans =

7.0108e-10

Your effort is highly appreciated!

Solution

  • By default, integral expects the function handle to take a vector input. In your code, the if-statement creates a complication since the condition will evaluate to true only if all elements of u satisfy it. So, if u is a vector that has elements both greater than and less than (1-lambda)/2, inv will never be assigned.

    There are two options:

    1. Put the if-statement in a for-loop and iterate over all of the elements of u.
    2. Use logical indexes for the assignment.

    The second option is faster for large element count and, in my opinion, cleaner:

    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;