Pourquoi les distributions t inverses pour les petites valeurs diffèrent-elles dans Matlab et R? - r, matlab, statistiques, distribution, probabilité

Je voudrais évaluer les inverses de Student t-fonction de distribution pour les petites valeurs, par exemple, 1e-18, dans Matlab. Le degré de liberté est de 2.

Malheureusement, Matlab revient NaN:

tinv(1e-18,2)
NaN

Cependant, si j'utilise la fonction intégrée de R:

qt(1e-18,2)
-707106781

Le résultat est raisonnable. Pourquoi Matlab ne peut-il pas évaluer la fonction pour cette petite valeur? Les résultats Matlab et R sont assez similaires à 1e-15, mais pour des valeurs plus petites, la différence est considérable:

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

Est-ce que quelqu'un sait quelle est la différence entre les algorithmes implémentés de Matlab et R, et si R donne des résultats corrects, comment pourrais-je calculer efficacement l'inverse t-distribution, dans Matlab, pour des valeurs plus petites?

Réponses:

5 pour la réponse № 1

Il semble que R "s qt peut utiliser un algorithme complètement différent que Matlab "s tinv. Je pense que vous et d'autres devriez signaler cette lacune à The MathWorks en déposant un demande de service. À propos, dans R2014b et R2015a, -Inf est retourné au lieu de NaN pour les petites valeurs (environ eps/8 et moins) du premier argument, p. C'est plus raisonnable, mais je pense qu'ils devraient faire mieux.

En attendant, il existe plusieurs solutions de contournement.

Cas spéciaux
Premièrement, dans le cas du Distribution t de Student, il y a plusieurs solutions analytiques simples à l'inverse CDF ou fonction quantile pour certains paramètres entiers de ν. Pour votre exemple de ν = 2:

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

qui retourne -7.071067811865475e+08. Au minimum, Matlab "s tinv devraient inclure ces cas particuliers (ils ne le font que pour ν = 1). Cela améliorerait probablement également la précision et la rapidité de ces solutions particulières.

Inverse numérique
le tinv fonction est basée sur la betaincinv une fonction. Il apparaît que c'est peut-être cette fonction qui est responsable de la perte de précision pour les petites valeurs du premier argument, p. Cependant, comme suggéré par l'OP, on peut utiliser la fonction CDF, tcdfet des méthodes de recherche de racines pour évaluer numériquement la CDF inverse. le tcdf la fonction est basée sur betainc, qui ne semble pas aussi sensible. fzero:

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

Cela retourne -7.071067811865468e+08. Notez que cette méthode n'est pas très robuste pour les valeurs de p proche de 1.

Solutions symboliques
Pour des cas plus généraux, vous pouvez profiter de maths symboliques et arithmétique à précision variable. Vous pouvez utiliser des identités en termes de Fonctions hypergéométriques gausiennes, ₂F₁, comme donné ici pour le CDF. Ainsi, en utilisant solve et 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)

le tinv fonction est basée sur la betaincinv une fonction. Il n'y a pas de fonction équivalente ou même une fonction bêta incomplète dans la boîte à outils Symbolic Math ou MuPAD, mais une fonction similaire ₂F₁ relation pour le fonction bêta incomplète peut être utilisé:

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)

Les deux schémas symboliques renvoient des résultats qui concordent avec -707106781.186547523340184 en utilisant la valeur par défaut de digits.

Je n'ai pas entièrement validé les deux méthodes symboliques ci-dessus, je ne peux donc pas garantir leur exactitude dans tous les cas. Le code doit également être vectorisé et sera plus lent qu'une solution entièrement numérique.

---