Fonction de répartition de la loi de Student

D'accord avec Gérard, vous ne trouverez pas de primitive.
De quoi avez-vous droit de vous servir pour faire le programme ?

"Beta incomplète"?. Dans mes classes, j'ai surtout des bétas complets.

(voir la réponse du Furet sur les STT, hier, dans "de l'intuition avant toute chose").

Bonsoir,

voici ma proposition.

1) Petit topo sur les fonctions Beta, Gamma et Beta incomplète ( n'en déplaise à RAJ .. :-) )


1.1. Fonction Gamma

La fonction Gamma est définie sous forme intégrale par

$$\Gamma(x)=\int_0^{+\infty} t^{x-1} e^{-t} dt$$

pour tout $x > 0$.

Pour $x=n$ entier naturel, on a

$$\Gamma(n)=(n-1)!$$

Pour $x=n+\frac{1}{2}$ où $n$ entier naturel

$$\Gamma(n+\frac{1}{2})=\frac{(2n-1) !! \sqrt{\pi}}{2^n}=\frac{(2n-1)! \sqrt{\pi}}{2^{2n-1} (n-1)!}$$

où $n!!$ désigne la double factorielle.

Enfin, une valeur particulière :

$$\Gamma(1/2)=\sqrt{\pi}$$


1.2. Fonction Beta

La fonction Beta est définie sous forme intégrale par

$$B(a,b)=\int_0^1 t^{a-1} (1-t)^{b-1} dt $$

pour $a >0$ et $b>0$. On peut montrer par un simple changement de variable que

$$B(a,b)=B(b,a)$$.

D'autre part, on a

$$B(a,b)= \frac{\Gamma(a) \Gamma(b)}{\Gamma(a+b)}$$


1.3. Fonction Beta incomplète

La fonction Beta incomplète est définie sous forme intégrale par

$$B(x;a,b)=\int_0^x t^{a-1} (1-t)^{b-1} dt$$

pour $0 \leq x \leq 1$, $a >0$ et $b>0$. En utilisant Chasle et un simple changement de variable, on peut montrer que

$$B(x;a,b)=B(a,b)-B(1-x;b,a)$$

On possède aussi le développement en série entière de la fonction Beta incomplète

$$B(x;a,b)=x^a \left( \frac{1}{a} + \frac{(1-b)}{(a+1)} x + ... + \frac{(1-b)(2-b)...(n-b)}{n! (a+n)} x^n + .... \right)$$

Ainsi, on peut avoit une bonne approximation de $B(x;a,b)$ en tronquant ce développement à un certain ordre. Il faut remarquer ( on s'en servira plus tard) que cette approximation est d'autant plus exacte si $x$ est petit.

2) Transformation du problème initial

On cherche à calculer

$$F(x)=\int_{-\infty}^x \frac{1}{\sqrt{d} B(1/2,d/2)} (1+\frac{t^2}{d})^{-\frac{d+1}{2}} dt$$

qui est la fonction de répartition d'une loi de Student à $d$ degrés de liberté.
Comme $F(+\infty)=1$, on a

$$\int_{-\infty}^{+\infty} \frac{1}{\sqrt{d} B(1/2,d/2)} (1+\frac{t^2}{d})^{-\frac{d+1}{2}} dt=1$$

et par parité de l'intégrande

$$\int_{-\infty}^0 \frac{1}{\sqrt{d} B(1/2,d/2)} (1+\frac{t^2}{d})^{-\frac{d+1}{2}} dt=\frac{1}{2}$$

Donc

$$F(x)=\int_{-\infty}^x \frac{1}{\sqrt{d} B(1/2,d/2)} (1+\frac{t^2}{d})^{-\frac{d+1}{2}} dt=\int_{-\infty}^0 \frac{1}{\sqrt{d} B(1/2,d/2)} (1+\frac{t^2}{d})^{-\frac{d+1}{2}} dt + \int_0^x \frac{1}{\sqrt{d} B(1/2,d/2)} (1+\frac{t^2}{d})^{-\frac{d+1}{2}} dt$$

Soit

$$F(x)=\frac{1}{2}+\int_0^x \frac{1}{\sqrt{d} B(1/2,d/2)} (1+\frac{t^2}{d})^{-\frac{d+1}{2}} dt$$

Il faut donc calculer

$$\int_0^x \frac{1}{\sqrt{d} B(1/2,d/2)} (1+\frac{t^2}{d})^{-\frac{d+1}{2}} dt = \frac{1}{\sqrt{d} B(1/2,d/2)} \int_0^x (1+\frac{t^2}{d})^{-\frac{d+1}{2}} dt = \frac{1}{\sqrt{d} B(1/2,d/2)} * I$$

Maintenant, on va transformer l'intégrale $I$ :

$$I=\int_0^x (1+\frac{t^2}{d})^{-\frac{d+1}{2}} dt =\frac{\sqrt{d}}{2} \int_0^{\frac{x^2/d}{1+x^2/d}} u^{-1/2} (1-u)^{d/2-1} dt$$

en faisant le changement de variable $\frac{t^2/d}{1+t^2/d}=u$ (à vous de le vérifier).

Donc, en utilisant le fonctions définies en 1), on a

$$I=\frac{\sqrt{d}}{2} B( \frac{x^2/d}{1+x^2/d}; 1/2, d/2)$$

En recollant tous les morceaux, on a ainsi

$$F(x)=\frac{1}{2} + \frac{1}{2} \frac{B(z;1/2,d/2)}{B(1/2,d/2)}$$

où $z=\frac{x^2/d}{1+x^2/d}$.

On a deux problèmes :

- calculer $B(1/2,d/2)$
- calculer $B(z;1/2,d/2)$


3) Calcul de $B(1/2,d/2)$

On utilise l'égalité de 1) :

$$B(1/2,d/2)=\frac{\Gamma(1/2) \Gamma(d/2)}{Gamma((d+1)/2)}$$

Alors en utilisant les résultat de 1) sur les fonctions $\Gamma$ pour $n$ et $n+\frac{1}{2}$ on trouve :

si $d=2n$ pair

$$B(1/2,d/2)=\frac{(n-1)! (n-1)!}{(2n-1)! 2^{2n-1}}$$

et si $d=2n+1$ impair

$$B(1/2,d/2)=\frac{\pi (2n-1)!}{2^{2n-1} (n-1)! n!}$$

Ce problème est donc réglé. Il faudrait que vos vérifiiez mes calculs ( je n'en suis pas sûr) en utilisant les résultats que j'ai cités au début et en utilisant les références que je vous donne à la fin.


4) Calcul de $B(z;1/2,d/2)$

C'est là qu'intervient la seule approximation (pour l'instant, tous les calculs sont exacts). On va utiliser le développement proposé en 2). Il est facile à coder, puisqu'on peut faire un programme itératif, qui à chaque itération rajoute un terme en plus du développement. Je ne détaille pas, vous êtes en école d'info !! :-)

Deux petites remarques cependant :

- il y a deux possibilités pour la partie itérative : soit vous fixez le nombre d'itération au départ ( par exemple 100, ce qui fait que vous prendrez le développement jusqu'à l'ordre 100), soit, et c'est la méthode que je vous conseille, vous arrêtez l'algorithme une fois que la différence entre la valeur obtenue à l'ordre $n$ et la valeur obtenue à l'ordre $n+1$ est suffisamment petite inférieure à $\epsilon$, personnellement je prendrais $\epsilon=10^{-l}$ où $l$ est le nombre de chiffres après la virgule de la valeur obtenue à l'ordre $n$.

- j'ai précisé plus haut que le développement était d'autant plus précis que $z$ était petit. Rappelons que $z$ est compris entre $0$ et $1$. On peut alors 'couper' la poire en deux, en utilisant la formule de passage de $B(z;a,b)$ à $B(1-z;b,a)$. En effet si $z$ est plus petit que $1/2$, on utilisera la formule de base. Sinon, si $z$ est plus grand que $1/2$, alors on utilisera

$$B(z;a,b)=B(a,b)-B(1-z;b,a)$$

et on pourra prendre le développement de $B(1-z;b,a)$ où $1-z$ sera plus petit que $1/2$. La limite de $1/2$ pour $z$ se traduit, en revenant à $x$, par une limite de $\sqrt{d}$ pour $x$.

En résumé (car je ne sais pas si j'ai été clair) :


- on calcule B=B(1/2,d/2) en utilisant les formules de 3)

- on a l'algo que je note BETA(x;a,b), algo itératif qui fait le développement en série, la règle d'arrêt des itérations étant définie comme ci-dessus ;

- si $x \qeq \sqrt{d}$ ( $z =\frac{x^2/d}{1+x^2/d} 1/2$) on fait $a=d/2$, $b=1/2$ et donc

$$F(x)=1-\frac{1}{2} \frac{BETA(1-z;a,b)}{B}$$

5) Une dernière remarque

Comme l'a signalé Gérard, si $d=1$, on tombe sur une loi de Cauchy dont on connaît la fonction de répartition :

$$F(x)=\frac{1}{2} + \frac{1}{\pi} \textrm{arctan}(x)$$

Donc faites tout d'abord un test pour savoir si $d=1$, et si c'est le cas, vous avez directement l'expression !

6) Quelques références électroniques


Propriétés des lois Gamma et Beta :

\lien{http://functions.wolfram.com/GammaBetaErf/}

\lien{http://mathworld.wolfram.com/IncompleteBetaFunction.html}

\lien{http://mathworld.wolfram.com/GammaFunction.html}


Voilà, a priori cette technique marche, il me semble que c'est même à peu près de cette façon que la fonction est implantée dans Matlab....

Amicalement,

ps : Ouf, j'ai fait le plus long poste de ma vie ... :-)

Tout à fait RAJ, j'ai oublié la plus jolie !
Cependant, ça fera un peu de boulot à ayobo de regarder ça, j'ai quand même maché pas mal le travail ...

ps : la longueur du post fait qu'il reste quelques coquilles de Latex ... un $x \geq \sqrt{d}$ et un $\Gamma$. Si jamais un modérateur passe dans le coin ...( ça fait deux fois aujourd'hui, promis j'arrête)

Amicalement,

$$B(1/2,d/2)=\frac{\Gamma(1/2) \Gamma(d/2)}{\Gamma((d+1)/2)}$$

Kuja : J'ai un doute ! La fonction densité est très lisse; Une méthode d'intégration genre Simpson ne donne-t-elle pas de très bons résultats ?

Bien évidemment, je trouve ton travail formidable et intéressant. Ce n'était qu'un doute. Et je n'irai pas vérifier concrétement quelle méthode est la meilleure. Par contre, tes arguments sont forts, et ta méthode doit être très efficace.
Amicalement,