 |
Introduction à la Recherche OpérationelleM2 Statistique Science des Données (SSD) |
Anatoli Juditsky Franck Iutzeler |
1- Problème de régression parcimonieuse et Dantzig Selector¶
On considère le modèle de régression normale $$ y=X\theta+\sigma\xi,\;\;\xi\sim \mathcal{N}(0, I_m),$$ où $X\in \mathbb{R}^{m\times n}$ et $y\in \mathbb{R}^m$ sont les observables, et $\theta\in \mathbb{R}^n$ est le paramètre inconnu.
L'estimateur de $\theta$ de ``Dantzig Selector'' (cf. Candes, E., Tao, T. (2007). The Dantzig selector: Statistical estimation when $p$ is much larger than $n$. The Annals of Statistics, 2313-2351) peut être utilisé pour estimer $\theta$ dans le cas d'un modèle surparamétré, quand la dimension $n$ de $\theta$ est supérieure a la dimension de l'observation $y$.
Dans ce cas l'estimateur ${\theta}_{DS}$ s'écrit comme une solution du probleme d'optimisation $$\min_{\theta\in \mathbb{R}^n} \|\theta\|_1,~~~\;\mbox{sous contrainte}\;\|X^T(X\theta-y)\|_\infty\leq \kappa\sigma,$$ où $\kappa>0$ est un hyper-paramètre.
La valeur de $\kappa$, préconisée dans la literature, est $c q_{\mathcal{N}}\left(1-{\alpha\over 2m}\right)$, où $\alpha\in (0,1)$ est le niveau de risque choisi (par ex. $\alpha=.05$) et $q_\mathcal{N}(p)$ est la $p$-quantile de la normale standardisée, et $c=\max_j\|[X]_j\|_2=\max_j\sqrt{[X]_j^T[X]_j}$ est la norme maximale de colonne de la matrice $X$.
Votre objectif dans cet exercice sera d'implementer l'estimateur ${\theta}_{DS}$ en utilisant CVXPY.
Question: Vérifier que le problème et les contraintes peuvent se formuler via des fonctions disponibles pour CVXPY.
2- Un petit exemple¶
import cvxpy as cp
import numpy as np
import pandas as pd
# Exemple jouet
n = 5
m = 3
sigma = 0.1
X = np.random.randn(m,n)
theta_true = np.array([1,0,0,0,5])
xi = np.random.randn(m)
y = np.dot(X,theta_true) + sigma*xi
X
y
theta_true
Question: Trouver l'estimateur $\theta_{DS}$ à partir de $X$ et $y$ par résolution du problème d'optimisation via CVXPY avec $\kappa$ fixé à 0.2.
theta = cp.Variable(n)
F = cp.norm1(theta)
objective = cp.Minimize(F)
kappa = 0.2
contrainte = cp.norm_inf( X.T@( X @ theta - y ) ) <= kappa*sigma
problem = cp.Problem(objective , constraints = [contrainte])
problem.solve(verbose=True,solver="ECOS")
problem.status
theta_ds= theta.value
theta_ds
print("True \t Dantzig selector")
for i in range(len(theta_true)):
print('{} \t {:6.3f}'.format(theta_true[i],theta_ds[i]))
3- Fonction "Dantzig Selector"¶
Question: Écrivez une fonction
DSelectqui fait appel aCVXPYpour calculer l'estimation ${\theta}_{DS}$.
Cette fonction doit sortir un tuple avec les elements
coef, vecteur des coefficients de regressionresid, vecteur $y-X{\theta}_{DS}$ de résidusstatus, le statut de sortie du solver
L'appel à cette fonction devra être:
DSelect(X, y, sigma = 1, c = 1, verb = False)
où
Xetysont les observablessigmaest une estimation de $\sigma$cest le paramètre réel tel que la valeur de $\kappa$ dans {DS} est $ \kappa=c\,q_{\mathcal{N}}\left(1-{\alpha\over 2m}\right).$
import scipy.stats
def DSelect(X, y, sigma = 1, c = 1, verb = False):
m,n = X.shape
alpha = 0.05
kappa = c*scipy.stats.norm.ppf(1-alpha/m)
theta = cp.Variable(n)
# Function
F = cp.norm1(theta)
objective = cp.Minimize(F)
# Constraints
contrainte = cp.norm_inf( X.T@( X @ theta - y ) ) <= kappa*sigma
problem = cp.Problem(objective , constraints = [contrainte])
problem.solve(verbose=verb,solver="ECOS")
theta_ds = theta.value
residual = y - np.dot(X,theta_ds)
sol_status = problem.status
return theta_ds,residual,sol_status
Test quand on connait le vrai theta¶
n = 5
m = 5
X = np.random.randn(m,n)
theta_true = np.random.randn(n)
sigma = 0.005
y = np.dot(X,theta_true) + sigma*np.random.randn(m)
Avec la valeur de sigma prise ci-dessus, vous devriez pouvoir appeler votre fonctions comme ci-dessous.
f = DSelect(X, y, sigma=sigma)
Question: Testez votre fonction sur cet exemple jouet. En particulier, comparez la valeur de
theta_dsà la vraie valeurtheta_true.
print("True \t Dantzig selector")
for i in range(len(theta_true)):
print('{:6.3f} \t {:6.3f}'.format(theta_true[i],f[0][i]))
Example du papier de Candes/Tao¶
import random
n = 256
m = 72
S = 8
S_set = random.sample(range(n),k=S)
X = np.random.randn(m,n)
theta_true = np.zeros(n)
theta_true[S_set] = np.random.randn(S)
sigma = 1/3.0*np.sqrt(S/m)
y = np.dot(X,theta_true) + sigma*np.random.randn(m)
f = DSelect(X, y, sigma= sigma)
Question: Testez votre fonction sur l'exemple du papier. Comparez la valeur de
theta_dsà la vraie valeurtheta_true. Essayer de changersigma.
f = DSelect(X, y, sigma= sigma)
np.linalg.norm(theta_true-f[0],1)
Support des solutions¶
La valeur cherchée theta_true est très parcimonieuse (elle contient beaucoup de zéros). Les seules valeurs non-null sont aux positions suivantes.
S_set
theta_true
Question: Comparez le nombre et l'emplacement des zéros dans
theta_dsà ceux danstheta_true. Essayez de changerc.
f[0]
import matplotlib.pyplot as plt
plt.plot(theta_true,'ro',label='Vrai theta')
plt.plot(f[0],'b*',label='Dantzig Selector')
plt.legend();
