为什么小值的逆t分布在Matlab和R中有所不同？

Question

我想在Matlab中为小值(例如1e-18)评估逆学生t-分布函数.自由度是2.

不幸的是,Matlab回归NaN:

tinv(1e-18,2)
NaN

但是,如果我使用R的内置函数:

qt(1e-18,2)
-707106781

结果是明智的.为什么Matlab不能评估这个小值的函数？Matlab和R的结果与1e-15非常相似,但对于较小的值,差异相当大:

tinv(1e-16,2)/qt(1e-16,2) = 1.05

有谁知道Matlab和R的实现算法有什么不同,如果R给出正确的结果,我怎么能在Matlab中有效地计算较小值的逆t-分布？

Answer 1

似乎R qt可能使用与Matlab 完全不同的算法tinv.我认为您和其他人应该通过提交服务请求向The MathWorks报告此缺陷.顺便说一句,在R2014b和R2015a,-Inf返回而不是NaN为小值(大约eps/8第一个参数的少), p.这更明智,但我认为他们应该做得更好.

在此期间,有几种解决方法.

特殊情况
首先,在学生t分布的情况下,对于ν的某些整数参数,对于反CDF或分位数函数有几种简单的解析解.对于你的例子ν = 2:

% for v = 2
p = 1e-18;
x = (2*p-1)./sqrt(2*p.*(1-p))

返回-7.071067811865475e+08.至少,Matlab tinv应该包括这些特殊情况(它们只对ν = 1 这样做).它也可能提高这些特定解决方案的准确性和速度.

数字反转
该tinv功能基于该betaincinv功能.似乎这个函数可能导致第一个参数的小值的精度损失,p.然而,如OP所建议的,可以使用CDF函数tcdf和根寻找方法以数字方式评估逆CDF.该tcdf功能基于betainc,似乎不敏感.使用fzero:

p = 1e-18;
v = 2
x = fzero(@(x)tcdf(x,v)-p, 0)

这回来了-7.071067811865468e+08.请注意,对于p接近的值,此方法不是非常健壮1.

符号解决方案
对于更一般的情况,您可以利用符号数学和变量精度算法.您可以使用高斯超几何函数(₂ F _1)中的标识,如此处给出的CDF.因此,使用solve和hypergeom:

% Supposedly valid for or x^2 < v, but appears to work for your example
p = sym('1e-18');
v = sym(2);
syms x
F = 0.5+x*gamma((v+1)/2)*hypergeom([0.5 (v+1)/2],1.5,-x^2/v)/(sqrt(sym('pi')*v)*gamma(v/2));
sol_x = solve(p==F,x);
vpa(sol_x)

该tinv功能基于该betaincinv功能.在符号数学工具箱或MuPAD中没有等效函数甚至不完整的Beta函数,但可以使用不完整Beta函数的类似₂ F ₁关系:

p = sym('1e-18');
v = sym(2);
syms x
a = v/2;
F = 1-x^a*hypergeom([a 0.5],a+1,x)/(a*beta(a,0.5));
sol_x = solve(2*abs(p-0.5)==F,x);
sol_x = sign(p-0.5).*sqrt(v.*(1-sol_x)./sol_x);
vpa(sol_x)

两个符号方案都返回同意-707106781.186547523340184使用默认值的结果digits.

我没有完全验证上面的两种符号方法,所以在所有情况下我都不能保证它们的正确性.代码也需要进行矢量化,并且比完全数值解决方案慢.