Décorrélation de variables aléatoires - correction

On construit des variables corrélées gaussiennes et on cherche à construire des variables décorrélées en utilisant le calcul matriciel. (correction)

Création d’un jeu de données

Q1

La première étape consiste à construire des variables aléatoires normales corrélées dans une matrice \(N \times 3\). On cherche à construire cette matrice au format numpy. Le programme suivant est un moyen de construire un tel ensemble à l’aide de combinaisons linéaires. Complétez les lignes contenant des .....

[1]:
import random
import numpy as np


def combinaison():
    x = random.gauss(0, 1)  # génère un nombre aléatoire
    y = random.gauss(0, 1)  # selon une loi normale
    z = random.gauss(0, 1)  # de moyenne null et de variance 1
    x2 = x
    y2 = 3 * x + y
    z2 = -2 * x + y + 0.2 * z
    return [x2, y2, z2]


li = [combinaison() for i in range(0, 100)]
mat = np.matrix(li)
mat[:5]
[1]:
matrix([[ 0.17179957, -0.80900898, -1.53083278],
        [-0.48698277, -0.95015277,  1.66842305],
        [ 0.53400158,  0.42468545, -2.41318128],
        [ 2.97269198,  9.62935172, -5.0672538 ],
        [ 0.88952227,  2.09758362, -2.36428305]])

Q2

[2]:
npm = mat
t = npm.transpose()
a = t @ npm
a /= npm.shape[0]
a
[2]:
matrix([[ 1.00968737,  3.07225686, -2.0111879 ],
        [ 3.07225686, 10.26264038, -5.16277675],
        [-2.0111879 , -5.16277675,  5.05453658]])

a est la matrice de covariance.

Corrélation de matrices

Q3

[3]:
cov = a
[4]:
var = np.array([cov[i, i] ** (-0.5) for i in range(cov.shape[0])])
var.resize((3, 1))
varvar = var @ var.transpose()
varvar
[4]:
array([[0.99040557, 0.31065402, 0.4426555 ],
       [0.31065402, 0.09744081, 0.13884485],
       [0.4426555 , 0.13884485, 0.19784207]])
[5]:
cor = np.multiply(cov, varvar)
cor
[5]:
matrix([[ 1.        ,  0.95440895, -0.89026339],
        [ 0.95440895,  1.        , -0.71682496],
        [-0.89026339, -0.71682496,  1.        ]])

Q4

[6]:
def correlation(npm):
    t = npm.transpose()
    a = t @ npm
    a /= npm.shape[0]
    var = np.array([cov[i, i] ** (-0.5) for i in range(cov.shape[0])])
    var.resize((3, 1))
    varvar = var @ var.transpose()
    return np.multiply(cov, varvar)


correlation(npm)
[6]:
matrix([[ 1.        ,  0.95440895, -0.89026339],
        [ 0.95440895,  1.        , -0.71682496],
        [-0.89026339, -0.71682496,  1.        ]])

Calcul de la racine carrée

Q6

Le module numpy propose une fonction qui retourne la matrice \(P\) et le vecteur des valeurs propres \(L\) :

L,P = np.linalg.eig(a)

Vérifier que \(P'P=I\). Est-ce rigoureusement égal à la matrice identité ?

[7]:
L, P = np.linalg.eig(a)
[8]:
P.transpose() @ P
[8]:
matrix([[ 1.00000000e+00,  2.91433544e-16, -5.55111512e-17],
        [ 2.91433544e-16,  1.00000000e+00, -2.77555756e-17],
        [-5.55111512e-17, -2.77555756e-17,  1.00000000e+00]])

C’est presque l’identité aux erreurs de calcul près.

Q7

np.diag(l) construit une matrice diagonale à partir d’un vecteur.

[9]:
np.diag(L)
[9]:
array([[1.44438859e+01, 0.00000000e+00, 0.00000000e+00],
       [0.00000000e+00, 1.60670907e-03, 0.00000000e+00],
       [0.00000000e+00, 0.00000000e+00, 1.88137173e+00]])

Q8

Ecrire une fonction qui calcule la racine carrée de la matrice \(\frac{1}{n}M'M\) (on rappelle que \(M\) est la matrice npm). Voir aussi Racine carrée d’une matrice.

[10]:
def square_root_matrix(M):
    L, P = np.linalg.eig(M)
    L = L**0.5
    root = P @ np.diag(L) @ P.transpose()
    return root


root = square_root_matrix(cov)
root
[10]:
matrix([[ 0.30467222,  0.77731521, -0.55914512],
        [ 0.77731521,  2.95469126, -0.96344225],
        [-0.55914512, -0.96344225,  1.9528626 ]])

On vérifie qu’on ne s’est pas trompé.

[11]:
root @ root - cov
[11]:
matrix([[ 1.77635684e-15,  1.77635684e-15, -1.33226763e-15],
        [ 1.33226763e-15, -1.06581410e-14,  4.44089210e-15],
        [-1.33226763e-15,  4.44089210e-15, -8.88178420e-16]])

Décorrélation

[12]:
np.linalg.inv(cov)
[12]:
matrix([[ 577.59780655, -117.8510365 ,  109.45002008],
        [-117.8510365 ,   24.24634204,  -22.12707857],
        [ 109.45002008,  -22.12707857,   21.14682284]])

Q9

Chaque ligne de la matrice \(M\) représente un vecteur de trois variables corrélées. La matrice de covariance est \(V=\frac{1}{n}M'M\). Calculer la matrice de covariance de la matrice \(N=M V^{-\frac{1}{2}}\) (mathématiquement).

Q10

Vérifier numériquement.

Simulation de variables corrélées

Q11

A partir du résultat précédent, proposer une méthode pour simuler un vecteur de variables corrélées selon une matrice de covariance \(V\) à partir d’un vecteur de lois normales indépendantes.

Q12

Proposer une fonction qui crée cet échantillon :

[13]:
def simultation(N, cov):
    # simule un échantillon de variables corrélées
    # N : nombre de variables
    # cov : matrice de covariance
    # ...
    return M

Q13

Vérifier que votre échantillon a une matrice de corrélations proche de celle choisie pour simuler l’échantillon.


Notebook on github