Problem 2: TINV

"tinv calls betainv. But for tinv(.5*(1+eps),10) for example, you're asking to find the quantile very close to 1 for a beta CDF that is very steep at 1. The current incarnation of betainv warns in this case"

p=0.5*(1+eps);
a=tinv(p,10)
er=tcdf(a,10)-p

Warning: BETAINV did not converge for a = 5, b = 0.5, p = 1.

a =

  9.5367e-007


er =

  3.7104e-007

There seems to be a problem with tcdf, too:

tcdf(1e-8,10)-0.5

ans =

     0

This number should be between

1e-8*[tpdf(0,10) tpdf(1e-8,10)]

ans =

  1.0e-008 *

    0.3891    0.3891

This is because tcdf calls betainc, and betainc is using the A&S formula that tcdf(x,v) = 0.5*betainc(v/(v+x^2),v/2,0.5);

Therefore, when x is less than v*sqrt(eps(1)), the number v+x^2 gets rounded to v, and the result is exactly 0.5. The problem is that x may not be as small as compared to v, whereas x^2 could be. This is the case when x=1e-8 and v=10.

Because of the same reason, the answer will be at most 8 digits accurate in the central region where the pdf is order 1. The situation worsens for higher values of v.

Solution: Write v/(v+x^2)=1-x^2/(v+x^2) and use AS 26.5.2 to see

tcdf(x,v)= 0.5*(1+sign(x).*betainc(x^2/(v+x^2),0.5,v/2));

Use this formula only in the central region as it would result in loss of precision in the tails.

x=1e-8; v=10;
0.5*(1+sign(x).*betainc(x^2/(v+x^2),0.5,v/2))-0.5-x*tpdf(x/2,v)

ans =

 -5.1151e-017

which is real tiny error, as we would expect from matlab. This change is now incorporated into tcdf and named mtcdf.m.

Next, we address tinv.m. Again, problems arise when we get close to the central region. On lines 51-54,

% q = p(k) - .5;
% df = v(k);
% z = betainv(1-2.*abs(q),df/2,0.5);
% x(k) = sign(q) .* sqrt(df ./ z - df);

when q is close to zero, z is close to 1, and also 1-z is close to zero. However, because of the way the problem is formulated, 1-z is not resolved well. The last line can be restated as

% x(k) = sign(q) .* sqrt(df .*(1-z)./z);

which explains the extra need for resolution of 1-z. Luckily, the problem can be fixed with the additional line

% oneminusz = betainv(2.*abs(q),0.5,df/2);

and the use of this variable in the calculation of x(k). These changes have been incorporated into tinv, and the renamed as mtinv.m.

The combination mtcdf with mtinv seem to be working well:

p=0.5+1e-8;

mtcdf(mtinv(p,10),10)-p

Error using ==> evalin
Undefined function or method 'mtinv' for input arguments of type 'double'.

as compared to

tcdf(tinv(p,10),10)-p

Just to check that tcdf is doing well, we run the loop:

x=linspace(-10,10,1000);
vlist=logspace(1,10,20);
for k=1:numel(vlist);
    v=vlist(k);
    p0=tcdf(x,v);
    p1=mtcdf(x,v);
    er(k)=max(abs(p1-p0));
end
max(abs(er))
max(abs(erx))

No difference detected when x is not too small, even when v is large.

Now, let us test tinv:

x=linspace(-10,10,1000);
vlist=logspace(1,10,20);
for k=1:numel(vlist);
    v=vlist(k);
    p1=mtcdf(x,v);
    imp=find(p1>1e-10 & p1 < 0.9999);
    p1=p1(imp); x=x(imp);
    x1=mtinv(p1,v);
    erx(k)=max(abs(x-x1));
end
max(abs(erx))

Errors may get big if very small p1 is tried to be inverted. (The reader may change p1 >1e-10 to p1>1e-20 to get erx = Inf). Some of this is because in tcdf, q=p-0.5 is defined and all the rest of the calculations are done with q. Instead, we try to use p directly whenever possible in mtinv. This impoves things to some extent, however mtinv still cannot invert mtinv(mtcdf(1e-18,1),1), for instance. It does invert mtinv(mtcdf(1e-15,1),1), although not very accurately.

Published with MATLAB® 7.4