Corrélation linéaire : génération données

Bonjour à tous

Je voudrais générer deux séries de données qui suivent une loi uniforme sur [0,1], corrélées linéaires.

Auriez-vous une méthode simple d'implémentation ?
On donne le coefficient de corrélation.

Merci
«1

Réponses

  • Je ne sais pas ce que tu appelles "corrélées linéaires" mais voir peut-être les copules archimédiennes classiques (que tu peux simuler dans R avec un package)
  • Si $X, Z$ sont indépendantes de loi $N(0,1)$ et si $Y=\cos \theta X+Z \sin \theta$ alors $\mathbb{E}(XY)=\cos \theta=r$ Soit $\Phi$ la fonction de répartion de la loi normale $N(0,1)$ et soit $U=\Phi(X)$ et $V=\Phi(Y)$. Soit un calul astucieux, soit un calcul laborieux montre que la corrélation entre $U$ et $V$ est $$\rho=\frac{6}{\pi}\arcsin \frac{r}{2}$$En résumé, pour simuler $(U,V)$ toutes deux uniformes sur (0,1) mais dépendantes avec corrélation $\rho$ tu dois prendre $r=\cos\theta=2\sin \frac{\pi}{6}\rho$ fabriquer $X,Z$ independantes $N(0,1)$ et écrire enfin
    $$U=\Phi(X),\ V=\Phi(X\cos \theta+Z\sin \theta)=\Phi( Xr+Z\sqrt{1-r^2}).$$
  • merci pour la méthode Gérard.
    bonne journée,
    Daniel S.
  • Bonjour

    Merci beaucoup pour votre réponse, M Letac.

    Et si vous pouviez détailler un peu "le calcul astucieux" qui montre la relation entre les deux coefficients de corrélation, je suis preneur.

    Cordialement
    Kezou
  • Autre question Gérard : comment calculer Psi(X), Psi n'étant pas explicitable ?

    Daniel
  • @danielsaada: De nos jours même Excel sait évaluer correctement $\Phi$, non ?
  • Steven, quel algorithme de calcul de Psi me proposes-tu
    pour implémenter la méthode de Gérard dans une sous-routine ?
  • Daniel,
    Pour ma part j'utiliserais un logiciel qui sait évaluer $\Phi$, quelle méthode utilise le logiciel j'en sais rien. Tu veux programmer dans quel langage ? Pourquoi ne pas utiliser une implémentation de $\Phi$ déjà faite ?
  • Steven,

    J'aime bien que mes programmes soient explicites et autonomes.
    Bonne soirée.
  • Daniel,

    Jusqu'à quel point ? Vas-tu implémenter la fonction exponentielle si tu en as besoin ?
  • ni l'exponentielle, ni l'addition !
  • Bonjour,

    Comme toujours, l'excellent Abramowitz et Stegun dispo en ligne fournit des approximations sympas, formules 26.2.16 à 26.2.19.

    Amicalement,
  • Alors tes progs ne sont pas autonomes, Daniel !
  • Je réponds avec retard à l'initiateur du fil qui me demandait la démonstration de $\frac{6}{\pi}\arg \sin\frac{r}{2}.$ Pardon, en particulier a Daniel Saada qui n'aime pas l'anglais, mais je n'ai pas le temps de traduire mon texte. Bah, les formules sont universelles.



    We start from the simplest idea: if $$\Phi(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^xe^{-\frac{u^2}{2}}du$$ and if $X\sim N(0,1)$ we have $\Phi (X)$ uniform on $(0,1)$. Let us denote \begin{equation}\label{TTT}T(x)=T(x)=2\sqrt{3}(\Phi(x)-1/2).\end{equation}
    \vspace{4mm}\noindent \vspace{4mm}\noindent \textbf{Proposition 4.1.} Let $(X,Y)$ be a centered Gaussian variable of $\mathbb{R}^2$ with covariance matrix $\left[\begin{array}{cc}1&r\\r&1\end{array}\right].$ Then
    \begin{eqnarray}\label{PHI}\mathbb{E}(\Phi(X)\Phi(Y))&=&\frac{1}{2}-\frac{1}{2\pi}\arg \cos\frac{r}{2}\\
    \label{THI}\mathbb{E}(T(X)T(Y))&=&\frac{6}{\pi}\arg \sin\frac{r}{2}\end{eqnarray}
    The proof of (\ref{PHI}) can be done by brute force and the computation of a four dimensional integral. We rather going to obtain Proposition 4.1 in a more interesting way after the following result.

    \vspace{4mm}\noindent \textbf{Theorem 4.2.} Let $(X,Y)$ be a centered Gaussian variable of $\mathbb{R}^2$ with covariance matrix $\left[\begin{array}{cc}1&r\\r&1\end{array}\right]$ and let $f$ be a real measurable function such that $\mathbb{E}_r(f(X))=0$ and $\mathbb{E}_r(f(X)^2)=1.$ Consider the Hermite polynomials $(H_n)_{n=0}^{\infty}$ defined by the generating function
    $$e^{xt-\frac{t^2}{2}}=\sum_{n=0}^{\infty}H_n(x)\frac{t^n}{n!}$$ and the expansion in orthogonal functions
    $$f(x)=\sum_{n=1}^{\infty}a_n\frac{H_n(x)}{\sqrt{n!}}.$$ Then for all $-1\leq r\leq 1$ we have $\sum_{n=1}^{\infty}a_n^2=1$ and
    \begin{equation}\label{AAA}\mathbb{E}(f(X)f(Y))=\sum_{n=1}^{\infty}a_n^2r^n.\end{equation}

    \vspace{4mm}\noindent \textbf{Proof.} Let us compute
    $$\mathbb{E}(e^{Xt-\frac{t^2}{2}}e^{Ys-\frac{s^2}{2}})= \sum_{n=0}^{\infty}\sum_{m=0}^{\infty}\frac{t^n}{n!}\frac{s^m}{m!}\mathbb{E}(H_n(X)H_m(Y))$$
    For this, write $r=\cos \alpha$ with $0\leq \alpha\leq \pi.$ If $X,Z$ are independent centered real Gaussian random variables with variance 1, then $Y=X\cos \alpha +Z\sin \alpha $ is centered with variance 1, $(X,Y)$ is Gaussian and $\mathbb{E}(XY)=\cos \alpha.$
    Therefore a simple calculation gives $$\mathbb{E}(e^{Xt-\frac{t^2}{2}}e^{Ys-\frac{s^2}{2}})=e^{ts\cos \alpha}$$
    This shows that
    $\mathbb{E}(H_n(X)H_m(Y))=0$ if $n\neq m$ and that $\mathbb{E}(H_n(X)H_n(Y))=n!\cos^n \alpha.$ From this we get the result. $\square$


    \vspace{4mm}\noindent \textbf{Corollary 4.3.} Let $p_n\geq 0$ such that $\sum_{n=1}^{\infty}p_n=1$ and consider the generating function $g(r)=\sum_{n=1}^{\infty}p_n r^n.$ Then $g(r)$ is the correlation the random variables
    $(f(X),f(Y))$ where $(X,Y)$ is centered Gaussian with covariance $\left[\begin{array}{cc}1&r\\r&1\end{array}\right]$ and where
    $$f(x)=\sum_{n=1}^{\infty}\epsilon _n\sqrt{p_n}\frac{H_n(x)}{\sqrt{n!}}$$
    with $\epsilon_n=\pm 1.$




    \vspace{4mm}\noindent \textbf{Proof of Proposition 4.1.} We apply Theorem 4.2 to the function $f=T$ defined by (\ref{TTT}). For this we have to compute $$\frac{a_n}{\sqrt{n!}}=\mathbb{E}(T(X)\frac{H_n(X)}{n!})$$
    Note that this is zero for even $n$ since $H_n$ and $T$ respectively even and odd functions. Thus we have to compute $p_{2n+1}> 0$ and $\epsilon_{2n+1}=\pm 1$ such that
    $$\epsilon_{2n+1}\frac{\sqrt{p_{2n+1}}}{\sqrt{(2n+1)!}}=\mathbb{E}(T(X)\frac{H_{2n+1}(X)}{(2n+1)!})$$
    To this purpose we watch the coefficient of $t^n$ in the power expansion of
    $$\mathbb{E}(T(X)e^{Xt-\frac{t^2}{2}})$$
    For this we need $$\mathbb{E}(\Phi(X)e^{Xt-\frac{t^2}{2}})=\Phi(\frac{t}{\sqrt{2}})=\frac{1}{2}+\frac{1}{2\sqrt{\pi}}\sum_{n=0}^{\infty}\frac{(-1)^n}{4^nn!}\frac{t^{2n+1}}{2n+1}$$$$\mathbb{E}(T(X)e^{Xt-\frac{t^2}{2}})=\sqrt{\frac{3}{\pi}}\sum_{n=0}^{\infty}\frac{(-1)^n}{4^nn!}\frac{t^{2n+1}}{2n+1}$$
    Therefore $$\epsilon_{2n+1}\frac{\sqrt{p_{2n+1}}}{\sqrt{(2n+1)!}}=\sqrt{\frac{3}{\pi}}\frac{(-1)^n}{4^nn!}\frac{1}{2n+1}$$
    which shows that $\epsilon_{2n+1}=(-1)^n.$
    To finish the proof we apply (\ref{AAA}) to $a_{2n+1}=\sqrt{p_{2n+1}}$ and we get
    $$\mathbb{E}(T(X)T(Y))=\frac{3}{\pi}\sum_{n=0}^{\infty}(\frac{1}{2})_n\frac{1}{4^nn!}\frac{r^{2n+1}}{2n+1}=\frac{6}{\pi}\arg \sin\frac{r}{2},$$ the last equality being easily checked. Of course (\ref{PHI}) is easily deduced from (\ref{THI}).$\square$
  • j'aime l'anglais mais je le comprends peu, nuance !
  • Une question pour Gérard : que signifie argsin et argcos ?

    Kuja, merci pour la référence !
  • Daniel, je voulais dire arcsinus et arccosinus. Millexcus.
  • mais Gérard, la dernière série de ton exposé en anglais
    n'est pas le développement en série de arcsin !

    Daniel
  • pas de $r\mapsto \arcsin r$, mais de $r\mapsto \frac{6}{\pi}\arcsin(r/2)$ certainement. Dérive les deux membres pour t'en convaincre ou bien développe $z\mapsto 1/\sqrt{1-z^2}$ en série entière, prends la primitive nulle en zéro et remplace $z$ par $r/2.$ Amicalement.
  • merci Gérard !
  • Gérard je suis désolé, ça ne marche pas,
    d'ailleurs j'ai soumis la série à un logiciel
    de calcul formel, il ne répond pas comme toi.
  • Est ce que tu interprètes correctement le symbole de Pochhhammer $(a)_n=a(a+1)(a+2)\ldots(a+n-1)?$

    Je détaillerai le calcul plus tard dans l'apres midi.
  • Gérard, ne te fatigue pas, j'avais cru lire 1/2^n
    et non le symbole de Pochhammer,
    Daniel
  • Bonsoir Gérard,

    J'épluche ligne à ligne la démonstration en anglais :
    peux-tu nous dire ce qu'est E indice r (ligne 11) ?

    Merci et bonne soirée,

    Daniel
  • Gérard, peux-tu me répondre, je suis au supplice !

    Daniel
  • C'est juste l'espérance, avec le $r$ en indice pour voir ça comme une fonction de la corrélation $r$.
  • merci de me délivrer Steven.
    bon dimanche,
    Daniel
  • Bonsoir à tous,

    pour ceux qui sont intéressés, j'ai rédigé la démonstration
    en anglais que Gérard Letac nous a transmise, et qui s'est
    révélée fort elliptique. Il reste deux points en suspens dans
    le 3) de la page 2 : je les ai marqués en gras.
    J'aimerais avoir votre avis.

    Cordialement,
    D.S
  • Bonsoir,

    Juste pour info, j'aurais tendance à penser après avoir écrit quelques lignes qu'il n'est pas nécessaire de passer par les polynômes d'Hermite, et que l'on peut simplement utiliser du calcul intégral "standard" pour calculer la corrélation.

    Amicalement,
  • je suis intéressé
  • Pour info, ce couplage de loi uniformes proposé par Gérard s'appelle copule de Gauss. Vous retrouverez ses formules dans ce doc : http://chess.uchicago.edu/ccehpe/hew_papers/Fall2007/Trivedi%2011_29_07.pdf
    Je ne crois pas qu'il y a les preuves, par contre elles sont probablement dans le livre de Nelsen.
  • Bonsoir,

    Calculs faits, ça peut aller très vite.
    La première étape est de démontrer très facilement que $E(G(X)G(Y)) = 12 E(F(X)F(Y)) - 3$
    Ensuite, la grosse astuce est de calculer la dérivée de $E(F(X)F(Y))$ par rapport à $\rho$ (la corrélation entre $X$ et $Y$).
    En utilisant le fait que
    \[\frac{\partial f_{X,Y}(x,y)}{\partial \rho}=\frac{\partial^2 f_{X,Y}(x,y)}{\partial x \partial y}\]
    puis Fubini, une IPP, Fubini et une IPP, on aboutit à
    \[\frac{\partial}{\partial \rho}E(F(X)F(Y)) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f_X(x)f_Y(y)f_{X,Y}(x,y)dxdy\]
    où je note $f_X$, $f_y$ et $f_{X,Y}$ les densités de $X$, $Y$ et $(X,Y)$ respectivement.
    Les calculs classiques sur le produit de lois normales donnent alors
    \[\frac{\partial}{\partial \rho}E(F(X)F(Y)) = \frac{1}{2\pi} \frac{1}{\sqrt{4-\rho^2}}\]
    On reconnaît alors la dérivée de $\arcsin(\rho /2)$, et en calculant la constante (pour $\rho=0$ par exemple), on retrouve le résultat.

    Amicalement,
  • merci beaucoup !
    Daniel
  • Bonjour,

    j'ai définitivement mis sur pied la démonstration de Gérard Letac
    (fichier joint) : c'est long mais instructif !
    Après un peu de repos, je m'attaque à celle de Kuja, que je remercie encore.
  • Il s'en passe des choses sur le forum dès qu'on tourne les talons...Ah, quel bosseur ce Daniel! Et bravo a Kuja qui raccourcit la méthode bourrin que j'avais mentionnée.
  • Ravi de te revoir Gérard !
    Ta démo est-elle de toi, sinon de qui ?

    Daniel
  • Bien sûr qu'elle est de moué: c'est plus amusant de bricoler soi même que de chercher dans la littérature. Gelfand dit que les grandes découvertes en mathématiques sont toujours faites par trois personnes en même temps...Ca m'étonnerait que ce que j'ai envoyé soit bien original. Amicalement
  • Gérard, ton nom a été mentionné ici à juste titre.

    Daniel
  • Comme je l'ai mentionné, ceci est la copule de Gauss, et j'aimerais bien savoir comment c'est généralisé au cas multivarié dans la littérature.
  • Par exemple, Stéphane, si tu cherches (U_1,U_2,U_3) uniformes et de matrice $R$ de corrélation donnée, il est faux de penser qu'il existe toujours une matrice $\Sigma$ de corrélation telle que si $(X_1,X_2,X_3)\sim N(0,\Sigma)$ alors $U_i=\Phi(X_i)$ fournit le résultat. Cherche sur Arkiv un article de 2010 de Luc Devroye et de ton serviteur, jamais publié ailleurs car toujours en préparation en vue d'un meilleur article de revue.
  • Ok mais ce que je me demandais c'est que vaut $R$ quand on donne $\Sigma$, est-ce ta preuve ou celle de Kuja qui s'adapte le mieux ?
    Enfin j'dis ça mais je n'y ai pas réfléchi un seul instant.
  • Salut Steven,

    Pour calculer $R$, il suffit de savoir calculer $R_{ij}$ pour tout couple $(i,j)$, donc on se ramène au cas $n=2$.
  • Par exemple, la matrice de corrélation $R$ de dimension 3 avec des -1/2 partout en dehors de la diagonale est inaccessible de cette facon gaussienne. C'est trop vague de dire qu'on se ramène au cas $n=2.$
  • Gérard : je parlais du "problème direct" à savoir calculer $R$ en fonction de $\Sigma$. Je suis bien conscient que le problème inverse est plus délicat.
  • Bonsoir à tous,

    maintenant qu'on sait construire deux va uniformes de coefficient
    de corrélation linéaire donné, quels problèmes intéressants peut-on
    poser ou se posent ?

    DS
  • Bonsoir,

    Je vais récupérer l'article de Gérard et Luc Devroye sur le champ, je ne savais pas que tu avais travaillé avec mon idole (à qui je fais souvent référence pour la simulation de variables aléatoires (entre autre) sur ce forum) !
    D'autre part, je plussoie Egoroff sur la question du problème direct. Mais pour revenir à mon "astuce", elle ne fonctionne que parce que les fonctions $\Phi$ (en reprenant la notation de Gérard) s'intègrent comme il faut. Dans le cas général, ce que je propose fait pschitt ...
    Enfin, pour relancer la question de Daniel, je pense qu'une idée pour se poser de nouveaux problèmes commencerait par aller regarder ce qui se fait du côté des copules avec la bilbio indiquée par Steven. Mais je n'y connais pas grand chose en copules 8-).

    Amicalement,
  • Une question classique avec les copules c'est le calcul de la corrélation de Kendall.
  • Je propose une solution alternative, qui me semble un peu plus simple que celles de Gérard et de Kuja, basée sur des manipulations de vecteur gaussiens. Ca revient à la "brute force and the computation of a four dimensional integral" mentionnée dans le texte de Gérard, mais en évitant au maximum de mettre les mains dans le cambouis.

    Je suppose $\rho > 0$, je note $\bar{\rho}=\sqrt{1-\rho^2}$. Toutes les lettres majuscules $X,Y,..$ représentent des gaussiennes standard indépendantes.

    On cherche à calculer $f(\rho)=\mathbb{E}(\Phi(X)\Phi(\rho X + \bar{\rho} Y))$. On écrit :
    \begin{align*}
    f(\rho) &= \mathbb{E}\left[ \mathbb{P} \left( Z \leq X, T \leq \rho X + \bar{\rho} Y | X,Y \right) \right] \\
    &= \mathbb{P} \left( Z \leq X, T \leq \rho X + \bar{\rho} Y \right) \\
    &= \mathbb{P} \left( X \geq Z, X \geq \rho^{-1} \left( T-\bar{\rho} Y \right) \right) \\
    &= \mathbb{P} \left( X \geq Z, X \geq cU \right)
    \end{align*}
    où $c^2 = \frac{1+\bar{\rho}^2}{\rho^2}$ est la variance de $\rho^{-1} \left( T-\bar{\rho} Y \right)$. Ensuite :
    \begin{align*}
    f(\rho) &= \mathbb{P} \left( X - Z \geq 0, X - cU \geq 0 \right) \\
    &= \mathbb{P} \left( \sqrt{2}V \geq 0, \sqrt{1+c^2} \left( \beta V + \bar{\beta} W \right) \geq 0 \right) \\
    &= \mathbb{P} \left( V \geq 0, \beta V + \bar{\beta} W \geq 0 \right)
    \end{align*}
    où $\beta=\left( 2(1+c^2) \right)^{-1/2}$ est la corrélation entre $X - Z$ et $X - cU$ et $\bar{\beta}=\sqrt{1-\beta^2}$.

    L'invariance par rotation de la loi de $(V,W)$ et un petit dessin montrent que la dernière probabilité vaut $\frac{1}{2}-\frac{1}{2\pi} \mathrm{arccos} \, \beta$, et un calcul rapide (le seul !!!) montre que $\beta = r/2$.
Connectez-vous ou Inscrivez-vous pour répondre.