.ipynb_checkpoints/Law_LAR-checkpoint.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "pycharm": {
     "is_executing": false,
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Law of the LAR's knots\n",
    "\n",
    "In this notebook, we present the law of the knots of the LAR. This result has been given in \n",
    "\n",
    "> *Multiple Testing and Variable Selection along Least Angle Regression's path*, [arXiv:1906.12072v3](https://arxiv.org/abs/1906.12072v3). \n",
    "\n",
    "## The recursive formulation of the LAR\n",
    "\n",
    "We begin with a new recursive algorithm, equivalent to LAR, that produces the sequence $(i_k,\\lambda_k)$ of indexes $i_k$ and knots $\\lambda_k$ of the LAR. This program has been implemnted under the name `lar_rec` given in the file `multiple_spacing_tests.py`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from lars.multiple_spacing_tests import lar_rec, generate_variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `generate_variables` gives an istance of a high-dimensional regression problem with *predictors_number* the number of predictors, *sample_size* the sample size, *sparsity=0* the sparsity of the target vector (each nonzero entry is a Rademacher random variable $\\{\\pm 1\\}$), $sigma=1$ the noise level. \n",
    "\n",
    "The matrix $Z$ is equal to $Z=X X^\\top$ where $X$ is a design matrix given by drawing i.i.d. Gaussian rows. The matrix $R$ is given by $R=[Z, -Z; -Z, Z]$. The target vector is *beta_0* and the response variable is $Y$.\n",
    "\n",
    "Here, $p=250$ (predictors number) and $n=50$ (sample size)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "[Z,R,sigma,X,noise_correlation,sparsity,beta_0,Y] = generate_variables(predictors_number=250, \n",
    "                                                                       sample_size=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute the LAR using the recursive formulation. See also the notebook `multiple_spacing_tests.ipynb` for a comparison with `SciKitLearn` version of LAR.\n",
    "\n",
    "It returns the knots *lambas* and the indexes *indexes*. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAboAAAEoCAYAAAAuWTTNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAArKUlEQVR4nO3deXyU5b338c9vJpN9gyQQyELCJiC7AVm07opbXY7Fpa6tBaune+tST/Xp6elpn57neHo8dalbXUtRscpR6oZaFQVMENl3BMIWCJIFsud6/shAIwYSIJl7ZvJ9v17zysw919zzyyWvfL3u+7rvy5xziIiIRCuf1wWIiIh0JQWdiIhENQWdiIhENQWdiIhENQWdiIhEtRivCzhemZmZrqCgwOsyREQkxEpKSnY757LaaxfxQVdQUEBxcbHXZYiISIiZ2aaOtNOhSxERiWoKOhERiWoKOhERiWoRf45ORCRcNTQ0UFpaSm1trdelRLT4+Hhyc3MJBALH9HkFnYhIFyktLSUlJYWCggLMzOtyIpJzjvLyckpLSyksLDymfejQpYhIF6mtrSUjI0MhdxzMjIyMjOMaFSvoRES6kELu+B1vHyroREQkqinogJr6JpqatS6fiEhnuvnmm1mxYoXXZSjo5izdztB7Xmfj7mqvSxER6VLOOZqbmzt1n01NTYd977HHHmPYsGGd+n3HotsHXZ+0eAA27t7vcSUiIp3v888/Z+jQodx6662MHTuWX/3qV4wbN46RI0dy7733Hmz39NNPM3LkSEaNGsV1110HwI033siLL754sE1ycjIA7733HmeccQbXXHMNI0aMYN++fVx44YWMGjWK4cOHM3PmTABOP/10iouLeeihh7j99tsP7ufJJ5/ke9/7HgDPPvss48ePZ/To0UyfPv2IwXmsuv3lBYWZSQB8vnufx5WISDT75f8uZ8W2yk7d57C+qdx78Ynttlu9ejV/+tOfuPTSS3nxxRdZuHAhzjm+/vWv8/7775ORkcGvf/1r5s2bR2ZmJnv27Gl3nwsXLmTZsmUUFhYya9Ys+vbty2uvvQZARUXFl9peccUVTJw4kd/97ncAzJw5k7vvvpuVK1cyc+ZM5s2bRyAQ4NZbb+W5557j+uuvP4beOLyQjejMLN7MFprZZ2a23Mx+2UYbM7P7zWydmS0xs7FdXVd6YizpiQE2livoRCQ69evXjwkTJvDmm2/y5ptvMmbMGMaOHcuqVatYu3Yt77zzDldccQWZmZkA9OzZs919jh8//uB1bSNGjODtt9/mjjvu4IMPPiAtLe1LbbOysujfvz/z58+nvLyc1atXM3nyZObOnUtJSQnjxo1j9OjRzJ07lw0bNnT67x/KEV0dcKZzrtrMAsCHZvY359z8Vm3OBwYFHycDDwV/dqnCzCSN6ESkS3Vk5NVVkpJajlw557jrrruYPn36l96///7725zCHxMTc/CcnnOO+vr6r+wTYPDgwZSUlDBnzhzuuusuzj33XO65554v7evKK6/k+eefZ8iQIVx22WWYGc45brjhBn7zm9902u/alpCN6FyLAzM+AsHHoVMdLwGeDradD6SbWZ+urq0wQ0EnItHvvPPO44knnqC6uuVP8datWykrK+Oss87i+eefp7y8HODgocuCggJKSkoAeOWVV2hoaGhzv9u2bSMxMZFrr72Wn/70pyxatOgrbS6//HJefvllZsyYwZVXXgnAWWedxYsvvkhZWdnB7920qUMr7xyVkJ6jMzM/UAIMBB5wzi04pEkOsKXV69Lgtu2H7GcaMA0gPz//uOsqyEzipU+3UlPfREKs/7j3JyISjs4991xWrlzJxIkTgZbJJc8++ywnnngid999N6eddhp+v58xY8bw5JNP8p3vfIdLLrmE8ePHc9ZZZ31pFNfa0qVL+dnPfobP5yMQCPDQQw99pU2PHj0YNmwYK1asYPz48QAMGzaMf/u3f+Pcc8+lubmZQCDAAw88QL9+/Tr19zbnQn/9mJmlA38FvuecW9Zq+2vAb5xzHwZfzwVud86VHG5fRUVF7ngXXp392Ta+P+NTXv/hqQzJTj2ufYmIHLBy5UqGDh3qdRlRoa2+NLMS51xRe5/15PIC59xe4D1gyiFvlQJ5rV7nAtu6up7CDM28FBGJVqGcdZkVHMlhZgnA2cCqQ5rNBq4Pzr6cAFQ457bTxQoyEwFdSyciEo1CeY6uD/BU8DydD3jeOfeqmd0C4Jx7GJgDXACsA/YDN4WisJT4AJnJcRrRiUinc87pxs7H6XhPsYUs6JxzS4AxbWx/uNVzB9wWqppaK8xM1LV0ItKp4uPjKS8v11I9x+HAenTx8fHHvI9uf2eUAwoyknhvzS6vyxCRKJKbm0tpaSm7dulvy/E4sML4sVLQBRVkJrGrpJTqukaS49QtInL8AoHAMa+KLZ2n29/U+QDd81JEJDop6IIKDlxioPN0IiJRRUEXdOASA43oRESii4IuKDE2huzUeF1LJyISZRR0rRRkJmqlcRGRKKOga6UwM4nPyzWiExGJJgq6Vgoyktizr56KmraXohARkcijoGulQJcYiIhEHQVdKwevpdMlBiIiUUNB10p+z0TMYKNGdCIiUUNB10p8wE/ftAQduhQRiSIKukMUZiZpRCciEkUUdIdouZZu33GvfyQiIuFBQXeIgowkKmsb+WK/LjEQEYkGCrpDHJh5qcOXIiLRQUF3CC3XIyISXRR0h8jrmYjfZ7qWTkQkSijoDhHw+8jtkcAGjehERKKCgq4NBRlJOnQpIhIlFHRtKMxsCTpdYiAiEvkUdG0oyEhkX30Tu6rrvC5FRESOk4KuDYVZyQB8rtXGRUQinoKuDYUZusRARCRaKOja0Dc9noDfNPNSRCQKKOjaEOP3kdczUSM6EZEooKA7jMKMJF00LiISBUIWdGaWZ2bvmtlKM1tuZj9oo83pZlZhZouDj3tCVd+hCjJbgq65WZcYiIhEspgQflcj8BPn3CIzSwFKzOwt59yKQ9p94Jy7KIR1takwM4nahmZ2VtXSJy3B63JEROQYhWxE55zb7pxbFHxeBawEckL1/UdLqxiIiEQHT87RmVkBMAZY0MbbE83sMzP7m5mdGNrK/qHg4CoGupZORCSShfLQJQBmlgzMAn7onKs85O1FQD/nXLWZXQC8DAxqYx/TgGkA+fn5XVJnn9R44mJ8bNxd3SX7FxGR0AjpiM7MArSE3HPOuZcOfd85V+mcqw4+nwMEzCyzjXaPOOeKnHNFWVlZXVKrz2f0y0hko0Z0IiIRLZSzLg14HFjpnLvvMG2yg+0ws/HB+spDVeOhCnSJgYhIxAvlocvJwHXAUjNbHNz2cyAfwDn3MHAF8F0zawRqgKuch0sIFGYl8d7qXTQ1O/w+86oMERE5DiELOufch8AR08I59wfgD6GpqH2FGUnUNzWzbW8NeT0TvS5HRESOge6McgQHZ17q8KWISMQK+azLSNI/GHR3zlrKmPx0hvVNZVifVE7sm0ZWSpzH1YmISEco6I6gV2o8/+fiYXy0vpzFW/by6pLtB9/LSoljWJ9Uvn1KIV8b3DUzP0VE5PiZh3M9OkVRUZErLi4OyXdV1DSwcnslK7ZVsmJ7JR+u3U2zc3x055nE+HUUWEQklMysxDlX1F47jeiOQlpCgAn9M5jQPwOAN5bvYPozJby3ehdnD+vtcXUiItIWDUOOw5lDepGZHMfM4i1elyIiIoehoDsOAb+PfxqbwzuryiirqvW6HBERaYOC7jhNHZdHU7NjVslWr0sREZE2KOiO04CsZMYV9OD54i1E+sQeEZFopKDrBFeOy2fj7n0s3LjH61JEROQQCrpOcMGIbJLjYjQpRUQkDCnoOkFibAwXj+rLnKXbqaxt8LocERFpRUHXSa4al0dtQzOzF2/zuhQREWlFQddJRuamMSQ7hed1+FJEJKwo6DqJmXHluDyWlFawYlul1+WIiEiQgq4TXTo6h1i/T6M6EZEwoqDrRD2SYjlveDZ//XQrtQ1NXpcjIiIo6DrdlUV5VNQ08MbyHV6XIiIiKOg63aQBGeT2SNDhSxGRMKGg62Q+nzG1KI9568rZXL7f63JERLo9BV0XuOKkXMzghRKN6kREvKag6wJ90xM4bXAWMxZuYWellu8REfGSgq6L/PicwdTUN3L1o/MpU9iJiHhGQddFRuam8+S3xrOjorYl7LQwq4iIJxR0XWhcQU+evGk82ytquebRBeyqqvO6JBGRbkdB18XGF/bkTzeOY+sXNVzz6Hx2VyvsRERCSUEXAif3z+BPN42jNBh25Qo7EZGQUdCFyIT+GTx+YxGb9+znm48tUNiJiISIgi6EJg3I5IkbxrFx9z6ueXQBby7foXtiioh0sZAFnZnlmdm7ZrbSzJab2Q/aaGNmdr+ZrTOzJWY2NlT1hcqkgZk8fsM4yqpqmfZMCWN/9Ra3PbeI2Z9to0qrk4uIdLqYEH5XI/AT59wiM0sBSszsLefcilZtzgcGBR8nAw8Ff0aVUwZlsvDus5m/oZy/LdvBm8t38trS7cT6fZwyKJMpw7M5f3g2KfEBr0sVEYl45pzz5ovNXgH+4Jx7q9W2PwLvOedmBF+vBk53zm0/3H6KiopccXFxl9fblZqaHYs2f8Hry3bw+rIdbN1bQ0pcDFefnM9Nkwvok5bgdYkiImHHzEqcc0XttvMi6MysAHgfGO6cq2y1/VXgt865D4Ov5wJ3OOeKD/n8NGAaQH5+/kmbNm0KVeldzjnH4i17eWLe58xZuh0Dvj6qLzef2p9hfVO9Lk9EJGx0NOhCeegSADNLBmYBP2wdcgfebuMjX0li59wjwCPQMqLr9CI9ZGaMye/B/+T34PbzTuCJeRuZ+ckWXvp0K6cOymTa1/pzysBMzNrqKhEROVRIZ12aWYCWkHvOOfdSG01KgbxWr3OBbaGoLRzl9Uzk3otP5OM7z+Jn553Aqh1VXPf4Qq5+dD579tV7XZ6ISEQI5axLAx4HVjrn7jtMs9nA9cHZlxOAiiOdn+su0hID3HbGQD684wx+dcmJLNq8l8senMf6XdVelyYiEvZCOaKbDFwHnGlmi4OPC8zsFjO7JdhmDrABWAc8CtwawvrCXlyMn+smFjDjOxOorm3ksgfm8dH63V6XJSIS1jybddlZomHW5bHYsmc/33ryEzbu3se/XzaCqePy2v+QiEgU6ehkFN0ZJULl9Uxk1q2TmDggg9tnLeG3f1tFc3Nk/0+LiEhXUNBFsNT4AE/cOI5rTs7n4b+v57Y/L6KmXrcUExFpTUEX4QJ+H7++dDj/cuFQXl++g6senU+lbiUmInKQgi4KmBk3n9qfh689ieVbK7jlmRLqGjWyExEBBV1UOe/EbH53xUg+Wl/OT57/TOfsRETw4M4o0rUuH5tLWVUdv/3bKnqlxPOLi4bqLioi0q0p6KLQ9K/1Z2dlLU/M20h2WhzTvjbA65JERDyjoItCZsYvLhxGWVUd/z5nFVkpcVw2JtfrskREPKGgi1I+n3Hf1FHsqa7nZy8sISMpjq8NzvK6LBGRkNNklCgWF+Pnj9efxKDeKdzybAlLSvd6XZKISMjpFmDdQFllLZc9+BE1DU2MK+gBgGEcmKNi1nK4c2h2ClOGZzOwV4qH1YqIdExYL7zamRR0HbN+VzV3vbSUiv0tF5M7HAf+0zugsamZz8v3AzAgK4kpw7OZcmIfhuekatamiIQlBZ0ctR0Vtby5YgevL9vBgo17aGp25KQncN6J2Vxzcj4DeyV7XaKIyEEKOjkue/bV8/bKnbyxbAcfrN1Nj6QA8+44kxi/TuuKSHjQ6gVyXHomxTK1KI/HbxzH/1wzhp2Vdfx9zS6vyxIROWoKOmnXmUN6kZEUywvFpV6XIiJy1BR00q6A38dlY3J4e+VOyqvrvC5HROSoKOikQ75RlEdjs+Plxdu8LkVE5Kgo6KRDTshOYVRuGi8UbyHSJzCJSPeioJMO+0ZRHqt2VLF0a4XXpYiIdJiCTjrs4lF9iYvxaVKKiEQUBZ10WFpCgCnDs3ll8VZqG7SCuYhEBgWdHJWpRXlU1jbyxvIdXpciItIhxxR0Zhbb2YVIZJjYP4Oc9AReLNHhSxGJDMc6ovvlgSdmNrmTapEI4PMZV5yUy4frdrN1b43X5YiItOtYg+4NM5tmZt8Ezu/MgiT8XXFSLs7BLI3qRCQCHHXQmdnjwNeBq4ATnXP/0ulVSVjL65nIpAEZvFCyheZmXVMnIuHtqIPOOfdt4G7gF0CZmf2x06uSsDe1KI8te2qYv7Hc61JERI7omA5dOudqnHPznHO/d85N78hnzOwJMyszs2WHef90M6sws8XBxz3HUpuExpTh2aTEx/CirqkTkTB3VEFnZnlmNsXMfmpmT5nZ0SwE9yQwpZ02HzjnRgcf/3o0tUloxQf8XDyqL3OWbaeytsHrckREDqvdoDOz6Wb2kZntBdYANwPJwGzgmo5+kXPufWDPMdYpYWhqUR61Dc28tmS716WIiBxWR0Z0dwE/Ak4CXgXigSecc7Occ2s6uZ6JZvaZmf3NzE48XKPgjM9iMyvetUuLgXplVG4ag3snM2PhZnZU1OpmzyISlqy9P05mNtw5t6zV6ynA/6XlUOR/O+eaO/xlZgXAq8654W28lwo0O+eqzeyC4L4HtbfPoqIiV1x8NEdQpTM9/fHn3PPKcqDlFmEn9E5hcHZyy8/eKZyQnUJ6ou4vICKdz8xKnHNF7bWLaa9B65ALvn7dzN4F/gWYB0w85iq/vN/KVs/nmNmDZpbpnNvdGfuXrnHdhH4M65PKyu2VrNpRxZqdVbyyeBtVtY0AxPiMR68v4owhvTyuVES6q3aDri3OuTrgF2b2TGcVYmbZwE7nnDOz8bQcVtXc9TBnZhQV9KSooOfBbc45dlbWsWpHJffOXs7v317D6SdkYWYeVioi3dUxBd0BR3OOzsxmAKcDmWZWCtwLBIL7eRi4AviumTUCNcBVTid9IpKZkZ0WT3ZaPDefWsMvXl7GJ59/wfjCnu1/WESkkx1X0B0N59zV7bz/B+APISpHQuSKsbnc9+ZqHnl/g4JORDyhZXqkSyXE+rluYgFzV+1k/a5qr8sRkW5IQSdd7vqJ/Qj4fTz+4UavSxGRbkhBJ10uMzmOfxqbw6ySUsqr67wuR0S6GQWdhMS3T+lPXWMzz8zf5HUpItLNKOgkJAb2SuasIb14+uNN1DY0eV2OiHQjCjoJmZtP7c+effW8tGir16WISDeioJOQmdC/JyNy0njsgw1asFVEQkZBJyFjZtx8aiEbdu9j7qoyr8sRkW5CQSchdcGIPuSkJ/DoBxu8LkVEugkFnYRUwO/jpskFLNy4h8+27PW6HBHpBhR0EnJXjssjJS5GozoRCQkFnYRcSnyAa07OZ87S7WzZs9/rckQkyinoxBM3Ti7AZ8b3//IpCzfu8bocEYliCjrxRJ+0BH5z+Qi27NnP1D9+zFWPfMxH63ejlZlEpLNZpP9hKSoqcsXFxV6XIceopr6JGQs38/Df11NWVcf4gp58/6xBTB6YoYVaReSIzKzEOVfUbjsFnYSD2oYmZn6yhYfeW8+OylrG5qfzz2cO5LTBvfD7FHgi8lUKOolIdY1NvFBcykPvrWfr3hpyeyRw9fh8phblkZUS53V5IhJGFHQS0eobm3lzxQ6em7+ZjzeUE+MzzhuezTdPzmdifx3WFBEFnUSRdWXVzFi4mRdLSqmoaaB/ZhI3TS7g2gn9FHgi3VhHg06zLiXsDeyVzC8uGsaCn5/Ff35jFCkJAX7xynL+vmaX16WJSARQ0EnEiA/4+aeTcnlh+kQyk2N5dv5mr0sSkQigoJOIExvj48pxebyzaidb99Z4XY6IhDkFnUSkq8fn44AZCzSqE5EjU9BJRMrtkciZJ/TiL59sob6x2etyRCSMKegkYl07sR+7q+t4c8UOr0sRkTCmoJOIddqgLPJ6JvDs/E1elyIiYUxBJxHL5zOuGd+P+Rv2sHZnldfliEiYUtBJRJtalEus38dzmpQiIoehoJOIlpEcxwUjsplVUsr++kavyxGRMBSyoDOzJ8yszMyWHeZ9M7P7zWydmS0xs7Ghqk0i27UT+lFV18jsxdu8LkVEwlAoR3RPAlOO8P75wKDgYxrwUAhqkihwUr8eDMlO4Zn5m7Rwq4h8RciCzjn3PrDnCE0uAZ52LeYD6WbWJzTVSSQzM745oR/Lt1WyeMter8sRkTATTufocoAtrV6XBrd9hZlNM7NiMyvetUs39hW4bEwOSbF+3f9SRL4inIKurfVW2jwO5Zx7xDlX5JwrysrK6uKyJBIkx8Vw2dgcXl2yjS/21XtdjoiEkXAKulIgr9XrXECzC6TDrp3Qj7rGZl4sKfW6FBEJI+EUdLOB64OzLycAFc657V4XJZFjSHYqRf168NyCTTQ3a1KKiLSICdUXmdkM4HQg08xKgXuBAIBz7mFgDnABsA7YD9wUqtokelw7oR8/nLmYac+UMLBXMjk9EsjtkUBejwRy0hNJiPV7XaKIhFjIgs45d3U77zvgthCVI1Hq/BHZvL82h0WbvuDva8poaPryyC4jKZa0hABm4PcZPmt5tDxv2RYb4yM2xk9cjI/YGB9x/uDPGB9xgZbt8cGfB7YlBPyML+xJ79R4j35zETmckAWdSCjExfi5b+poAJqbHWVVdWzdu5/SL2qCj/1U1zXR3OxoanY0uwMPaApuq29spmJ/PXWNzdQ3NVPf2PKoa2ymrrGJusZm2rpcz+8zzhnam29OyGfygEx8vrbmV4lIqCnoJGr5fEZ2WjzZafGc1K/z9uuco74pGHwNzdQ2NFFR08D/fraNF0pKeX35DvplJHLN+HyuOCmXjOS4zvtyETlqFul3kigqKnLFxcVelyECQF1jE68v28FzCzazcOMeYv0+pgzP5qKRfRidl04vHdoU6TRmVuKcK2q3nYJOpGus3VnFcws289KiUiprW244nZ0az4jcNEblpjEiN52ROWn0SIr1uFKRyKSgEwkTtQ1NLNtawZLSCpaU7mVJaQUbdu87+P6Q7BR+fM5gzhnWGzOd1xPpqI4Gnc7RiXSx+ICfooKeFBX0PLitsraBZaUVfFZawYslW5j2TAmTBmTwLxcOY1jfVA+rFYk+GtGJeKyhqZk/L9jMf729hoqaBqaelMdPzhtMrxSdzxM5ko6O6MLpzigi3VLA7+OGSQX8/adn8K3Jhbz0aSln/Md7PPDuOmobmrwuTyTiKehEwkRaYoBfXDSMN390GpMGZvIfb6zmrP/8O59p6SGR46KgEwkzhZlJPHp9EX/+zsn4fHDVI/N5d3WZ12WJRCwFnUiYmjQgk1nfnUT/rCRufqqYF4q3tP8hEfkKBZ1IGOuVEs/M6ROZNCCDn724hAfeXUekTyATCTUFnUiYS46L4fEbxnHp6L78xxuruXf2cpq0DJFIh+k6OpEIEBvj476po+mVGs8j729gV1Ud/3XlaOIDWnZIpD0KOpEI4fMZP79gKL1T4/nVqysor17Io9cXkZYY8Lo0kbCmQ5ciEebbpxTyP1ePYfGWvVxw/weUbNrjdUkiYU1BJxKBLh7Vl5nTJ+DzwdQ/zuf+uWt13k7kMBR0IhFqTH4P5nz/VC4e2Yf73lrD1Y/MZ+veGq/LEgk7CjqRCJYSH+D3V43hvqmjWL6tgvN//z5/W7rd67JEwoqCTiQKXD42l9e+fyqFmUl897lF3PXSEvbXN3pdlkhYUNCJRImCzCReuGUS00/rz4yFW7jo/g9ZuFETVUQUdCJRJDbGx13nD+W5m0+mvqmZqX/8mLv/upSq2gavSxPxjIJOJApNHpjJmz/6Gt+aXMifF27m3P96n7krd3pdlognFHQiUSoxNoZ7Lh7GS9+dREp8DN9+qpjvz/iU8uo6r0sTCSmtMC7SDdQ3NvPQe+v5w7trSY6L4UfnDCa/ZyJJcTEkxvpJio0hMa7lZ0LAj89nXpcs0q6OrjCuoBPpRtbsrOKOWUv4dPPew7YJ+I2LR/Xl1tMHMrBXcuiKEzlKCjoRaVNzs2P9rmqq6hrZX9dEdV0j++sb2VffxP66Rjbv2c+sRaXUNTZzwYg+3Hb6QIb1TfW6bJGv6GjQ6abOIt2Mz2cM6p1yxDY/Pmcwj3+4kac/3sRrS7Zz9tBe/POZgxidlx6aIkU6UUgno5jZFDNbbWbrzOzONt4/3cwqzGxx8HFPKOsTkRYZyXHcPmUI8+44kx+fM5jiTV9w6QPzuO7xBcxZup2K/bpcQSJHyA5dmpkfWAOcA5QCnwBXO+dWtGpzOvBT59xFHd2vDl2KdL3qukaem7+JRz/YwO7qenwGI3PT+dqgTE4dnMXovHQCfk3iltAKx0OX44F1zrkNAGb2F+ASYMURPyUinkuOi2H6aQP41imFfLZlL++v3c2Ha3fxh3fXcf8760iOi2FC/wwmDshgVG4aw/qmkhirMyMSHkL5LzEH2NLqdSlwchvtJprZZ8A2WkZ3yw9tYGbTgGkA+fn5XVCqiLQl4PdRVNCTooKe/PicwVTUNPDx+t28v3Y3H6zdxdvBi9J9BoN7pzAiJ42RuWmMyE1nSHaKVkQXT4Qy6Nq6MOfQ46aLgH7OuWozuwB4GRj0lQ859wjwCLQcuuzkOkWkg9ISAkwZ3ocpw/sAsLOyliWlFSwt3cuSrRXMXVXGCyWlB9vHxfhIioshKXjNXmKsn6S4GFLiY7hsTC5nD+2Fma7hk84VyqArBfJavc6lZdR2kHOustXzOWb2oJllOud2h6hGETkOvVPjOWdYPOcM6w2Ac45tFbUsLd3Lmp3V7KtrZF99y2UN++ob2V/fxL66RtburGbO0h2ML+jJnRcMYWx+D49/E4kmoQy6T4BBZlYIbAWuAq5p3cDMsoGdzjlnZuNpmRVaHsIaRaQTmRk56QnkpCcwZfjh2zU0NTPzky38/u21XP7gR5w/PJufnXcC/bN0wbocv5AFnXOu0cz+GXgD8ANPOOeWm9ktwfcfBq4AvmtmjUANcJWL9CvaRaRdAb+Payf047IxOTz6wQYeeX8Db63YydXj8/n+WYPISonzukSJYLozioiEnV1Vdfz33DXMWLiF+Bgf3zqlkG9NLqRHUqzXpUkY0S3ARCTibdhVzf97czVzlu4gMdbPtRP6cfOphfRKife6NAkDCjoRiRprdlbx4LvrmP3ZNmL8Pq4al8f00waQk57gdWniIQWdiESdz3fv4+G/r2fWolKcg8vG5HDrGQMpzEzyujTxgIJORKLWtr01PPL+BmYs3ExTs+PGSQV8/+xBpMYHvC5NQkhBJyJRb1dVHfe9tZq/fLKFjKRYbj9vCFeclKuFY7uJjgad7sIqIhErKyWO31w+ktm3nUK/jCRun7WEyx6cx6LNX3hdmoQRBZ2IRLwRuWm8eMtEfn/laHZU1nL5gx/x4+cXU1ZZ63VpEgZ0e3ERiQpmxqVjcjhnWG8eeHcdj32wkVc/287g7GQG90phcHYKg3snM7h3CjnpCbqnZjeic3QiEpU2le/j2fmbWLWjirU7q9nRanSXFOtnYO8U0hICxPqNgN938BEbY8T6fWQkxzEgK5kBvZIoyEjSygthKBzXoxMRCZl+GUncfeGwg68r9jewtqyKNTurWbOzinVl1VTWNNDQ1Bx8OOobW57XNzWzt9Uq6maQ1yORgb2SGZCVRG6PROJifMTG+IiL8RN78LmP+ICf3qlx9EqJx69JMWFBQSci3UJaYuDgWnodUVPfxIbd1azftY91ZdWs31XN+rJqPly3m/rG5nY/7/cZvVLi6JMWT5/0BPqkxpOdFk9C7OFHhrF+H2kJAVITAqS1eiTG+nWo9Tgo6ERE2pAQ6+fEvmmc2DftS9ubmh1f7K+nvrGZ+sZm6oI/65uaqGtspqa+iZ2VdWyvqGHb3lp2VNawclslc1fupLah/YBsS4zP6J+VxH9fNYahfVI749frVhR0IiJHwe8zMpOPfjUF5xwVNQ1HHA3WNTZTUdNw8FHZ6vmsRaX800Mfcf9VYzg7uN6fdIyCTkQkBMyM9MT2V1/IO8z2GyYVcPNTxXznmWJ+fv5Qbj61UIczO0jX0YmIRIDeqfE8P30i5w/P5tdzVnLnrKUdOlcoCjoRkYiREOvnD1eP5XtnDmRm8Raue3wBX+yr97qssKegExGJID6f8ZNzT+D3V47m0817ufTBeawrq/a6rLCmoBMRiUCXjslhxrQJ7Ktr5LIH5/Ffb61he0WN12WFJQWdiEiEOqlfD16+bTIn9evB/e+sZfJv3+Hmp4p5d3UZTc2RfderzqRZlyIiESy3RyJP3jSezeX7mfHJZl4o3sLbK3eSk57ANSfn842iXHqlxHtdpqd0r0sRkShS39jMWyt28tyCTXy0vpwYn3HeidncMKmAcQU9ouqSBN3rUkSkG4qN8XHhyD5cOLIPG3ZVM2PhZmZ+soXXlm5naJ9UbpzUj0tG53Srm1RrRCciEuVq6pt4efFWnvroc1btqCI9McCV4/K49uR+5PVM9Lq8Y9bREZ2CTkSkm3DOsXDjHp76+HPeWL4T5xzDc9JIT4wlNT6G1IQAqfEBUhNiSI0PkJkcR1FBj2O65Vko6NCliIh8iZlxcv8MTu6fwba9Nfx5wWY+K91LRU0DpXv2U1nbcl/NhqYvD4CGZKcwoX8Gkwa0fDYtIeDRb3BsNKITEZGDnHPUNTZTWdPAli9qmL+hnI/Xl1O8aQ+1Dc34DIbnpDGxfwYjctMYkp1KYWaSJ2vv6dCliIh0mrrGJhZv3stH68v5eEM5n27+4uDILy7GxwnZKQzNTmVonxSG9Emlf2YSmclx+LowABV0IiLSZeoam1hXVs3K7VWs3F7Jqh2VrNxexZ5W996MjfHRNy2enB4J5KQnkJOeSE6PBPJ7JjK+sGML4B5JWJ6jM7MpwH8DfuAx59xvD3nfgu9fAOwHbnTOLQpljSIi0r64mK8uTOuco6yqjpXbK9m8Zz9bv6ihdG8NW7+o4d3Vu9hVVQdATnoC8+48M2S1hizozMwPPACcA5QCn5jZbOfcilbNzgcGBR8nAw8Ff4qISJgzM3qnxtM7te07sdQ2NLGjopaq2saQ1hXKEd14YJ1zbgOAmf0FuARoHXSXAE+7luOp880s3cz6OOe2h7BOERHpAvEBPwWZSSH/3lDe1DkH2NLqdWlw29G2ERER6bBQBl1bU28OnQnTkTaY2TQzKzaz4l27dnVKcSIiEp1CGXSlQF6r17nAtmNog3PuEedckXOuKCsrq9MLFRGR6BHKoPsEGGRmhWYWC1wFzD6kzWzgemsxAajQ+TkRETkeIZuM4pxrNLN/Bt6g5fKCJ5xzy83sluD7DwNzaLm0YB0tlxfcFKr6REQkOoX0Ojrn3Bxawqz1todbPXfAbaGsSUREolsoD12KiIiEnIJORESimoJORESiWsTf1NnMdgGbjnM3mcDuTignmqmP2qc+OjL1T/vUR+1r3Uf9nHPtXmMW8UHXGcysuCN3wO7O1EftUx8dmfqnfeqj9h1LH+nQpYiIRDUFnYiIRDUFXYtHvC4gAqiP2qc+OjL1T/vUR+076j7SOToREYlqGtGJiEhUU9CJiEhU6/ZBZ2ZTzGy1ma0zszu9riccmNkTZlZmZstabetpZm+Z2drgzx5e1uglM8szs3fNbKWZLTezHwS3q4+CzCzezBaa2WfBPvplcLv6qBUz85vZp2b2avC1+qcVM/vczJaa2WIzKw5uO+o+6tZBZ2Z+4AHgfGAYcLWZDfO2qrDwJDDlkG13AnOdc4OAucHX3VUj8BPn3FBgAnBb8N+N+ugf6oAznXOjgNHAlODSW+qjL/sBsLLVa/XPV53hnBvd6tq5o+6jbh10wHhgnXNug3OuHvgLcInHNXnOOfc+sOeQzZcATwWfPwVcGsqawolzbrtzblHweRUtf6hyUB8d5FpUB18Ggg+H+uggM8sFLgQea7VZ/dO+o+6j7h50OcCWVq9Lg9vkq3ofWAQ3+LOXx/WEBTMrAMYAC1AffUnwsNxioAx4yzmnPvqy3wO3A82ttql/vswBb5pZiZlNC2476j4K6Xp0Ycja2KbrLaRDzCwZmAX80DlXadbWP6fuyznXBIw2s3Tgr2Y23OOSwoaZXQSUOedKzOx0j8sJZ5Odc9vMrBfwlpmtOpaddPcRXSmQ1+p1LrDNo1rC3U4z6wMQ/FnmcT2eMrMALSH3nHPupeBm9VEbnHN7gfdoOe+rPmoxGfi6mX1OyymTM83sWdQ/X+Kc2xb8WQb8lZbTTUfdR9096D4BBplZoZnFAlcBsz2uKVzNBm4IPr8BeMXDWjxlLUO3x4GVzrn7Wr2lPgoys6zgSA4zSwDOBlahPgLAOXeXcy7XOVdAy9+dd5xz16L+OcjMksws5cBz4FxgGcfQR93+zihmdgEtx8r9wBPOuV97W5H3zGwGcDoty2HsBO4FXgaeB/KBzcA3nHOHTljpFszsFOADYCn/OL/yc1rO06mPADMbSctEAT8t/0P9vHPuX80sA/XRlwQPXf7UOXeR+ucfzKw/LaM4aDnN9mfn3K+PpY+6fdCJiEh06+6HLkVEJMop6EREJKop6EREJKop6EREJKop6EREJKop6EREJKop6EREJKop6EQikJmdbWbPeF2HSCRQ0IlEplHAp14XIRIJFHQikWkU8KmZxZnZk2b276blE0Ta1N2X6RGJVKNouWv7G8BjzrlnPa5HJGzpXpciESa4RNBuYBMw3Tn3sccliYQ1HboUiTzDaFliqhFo8rgWkbCnoBOJPKOAj2hZx+xPZtbb43pEwpqCTiTyjAKWOefWAHcAzwcPZ4pIG3SOTkREoppGdCIiEtUUdCIiEtUUdCIiEtUUdCIiEtUUdCIiEtUUdCIiEtUUdCIiEtX+Pzl5NZ2yAJFpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "[lambdas, indexes, correls, R, Z] = lar_rec(X,Y)\n",
    "\n",
    "k_max = np.size(indexes)\n",
    "\n",
    "plt.plot(range(k_max-1),lambdas[:k_max-1])\n",
    "plt.legend([r'recursive'])\n",
    "plt.tight_layout()\n",
    "plt.xlabel(r'$k$')\n",
    "plt.ylabel(r'$\\lambda_k$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The joint law of the knots\n",
    "\n",
    "Given the joint law of the knots proved in \n",
    "\n",
    "> *Multiple Testing and Variable Selection along Least Angle Regression's path*, [arXiv:1906.12072v3](https://arxiv.org/abs/1906.12072v3),\n",
    "\n",
    "one can compute the conditional law of $\\lambda_b\\|\\lambda_a,\\lambda_c$ whose cdf is denoted by $\\mathbb F_{abc}$. In particuler, our theoretical analysis shows that $\\hat\\alpha_{abc}=\\mathbb F_{abc}(\\lambda_b)$ is uniformly distributed under the null (i.e., the target vector is zero). \n",
    "\n",
    "One can compute $\\hat\\alpha_{abc}$ using the function `observed_significance_CBC`. **CBC** means *Component-By-Component* as in \n",
    "\n",
    "> *Fast algorithms for component-by-component construction of rank-1 lattice rules in shift-invariant reproducing kernel Hilbert spaces*, by D. Nuyens, R. Cools in *Mathematics of Computation*, 2006/\n",
    "\n",
    "It is a lattice rule to compute the conditional marginal distribution of $\\lambda_b$. \n",
    "\n",
    "Here, we compute $\\hat\\alpha_{136}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Observed significance is 0.5749611629559863\n"
     ]
    }
   ],
   "source": [
    "from lars.multiple_spacing_tests import observed_significance_CBC\n",
    "\n",
    "# compute the LAR\n",
    "lars = [lambdas, indexes, correls]\n",
    "\n",
    "# a (start), b (middle) and c (end)\n",
    "start = 1\n",
    "middle = 3\n",
    "end = 6\n",
    "\n",
    "# hat alpha \n",
    "alpha, __, __ = observed_significance_CBC(lars, sigma=sigma, start=start, end=end, middle=middle)\n",
    "\n",
    "print(\"Observed significance is %s\" % alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By Monte Carlo, one can observe the law of $\\hat\\alpha_{abc}$. This is done by the command \n",
    "\n",
    "> `empirical_law_first_knots(predictors_number=250,\n",
    "                              sample_size=50,\n",
    "                              end=5,\n",
    "                              start=1,\n",
    "                              iterations=5000,\n",
    "                              display=1,\n",
    "                              eval_points=99991)`\n",
    "                              \n",
    "over `iterations=5000` Monte Carlo simulations and `eval_points=99991` points of Lattice rule. It took us several hours (around $20h$) to run this command and we give you the result here. Of course, you can check that you get the same result by running the above command. Thank you for your time!\n",
    "\n",
    "![Knots Law](img/laws_final.png)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": []
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}