Law_LAR.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 \"*Multiple Testing and Variable Selection along Least Angle Regression's path*\" (Azaïs and De Castro, 2019). \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 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, 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": "iVBORw0KGgoAAAANSUhEUgAAAboAAAEoCAYAAAAuWTTNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXzU1b3/8ddnJssEskAWkkACAUQh7BAQxKUuuC8titrdqrW2VWsftr/azbZ2u723V3t7bWtt69ZaxWpvpS6torhUEQjIvpVVQoCERLIA2WbO748ZMMZAEkjmO5m8n4/HPGY7851Pvg+SN+f7/Z5zzDmHiIhIvPJ5XYCIiEhPUtCJiEhcU9CJiEhcU9CJiEhcU9CJiEhcS/C6gBOVnZ3tioqKvC5DRESibNmyZfucczkdtev1QVdUVERpaanXZYiISJSZ2Y7OtNOhSxERiWsKOhERiWsKOhERiWu9/hydiEisam5upqysjIaGBq9L6dUCgQAFBQUkJiYe1+cVdCIiPaSsrIy0tDSKioowM6/L6ZWcc1RVVVFWVsbw4cOPaxs6dCki0kMaGhrIyspSyJ0AMyMrK+uEesUKOhGRHqSQO3Enug8VdCIiEtcUdEBDc5BgSOvyiYh0p/vvv59HH33U6zJ0McoLq3fzxceW8+JXz+Tk3DSvyxER6THOOZxz+Hzd28cJBoP4/f4PvX7zzTd36/ccrz7fo8tKTQZgT40u/xWR+LN9+3bGjBnDl770JaZMmcIf//hHZs6cyZQpU5g7dy719fUALF26lNNOO42JEycyffp06urqePjhh7nllluObOvSSy/l1VdfBSA1NZW77rqLU089lUWLFnHnnXdSXFzMhAkT+NrXvgbA97//fX7+85+zfv16pk+f/oGaJkyYAMCyZcs466yzmDp1KhdccAG7d+/u9n3Q53t0+RkBQEEnIj3rB39fy7ry2m7dZvHgdL532dgO223cuJGHHnqIu+++mzlz5rBgwQL69+/Pz372M+655x7uvPNOrrnmGubNm8e0adOora0lJSXlmNs8cOAA48aN4+6776a6upobbriBDRs2YGbs37//A23HjBlDU1MTW7duZcSIEcybN4+rr76a5uZmbr31Vp555hlycnKYN28e3/72t3nwwQdPaL+01eeDblB6pEdXq6ATkfg0bNgwZsyYwbPPPsu6deuYNWsWAE1NTcycOZONGzeSn5/PtGnTAEhPT+9wm36/nyuvvPJI+0AgwI033sgll1zCpZde+qH2V199NU8++SR33nkn8+bNY968eWzcuJE1a9Ywe/ZsIHwIND8/v7t+7CP6fNAlJ/jJ7J+koBORHtWZnldP6d+/PxA+Rzd79mwef/zxD7y/atWqdi/hT0hIIBQKHXneeixbIBA4cl4uISGBJUuW8PLLL/PEE09w33338corr3xgW9dccw1z585lzpw5mBmjRo1i9erVjB07lkWLFnXbz9qePn+ODiA3PaBDlyIS92bMmMGbb77J5s2bATh48CCbNm1i9OjRlJeXs3TpUgDq6upoaWmhqKiIFStWEAqF2LlzJ0uWLGl3u/X19dTU1HDxxRfzi1/8ghUrVnyozciRI/H7/fzwhz/kmmuuAeCUU06hsrLySNA1Nzezdu3abv+5+3yPDsLn6RR0IhLvcnJyePjhh/n4xz9OY2MjAD/60Y84+eSTmTdvHrfeeiuHDh0iJSWFBQsWMGvWLIYPH8748eMZN24cU6ZMaXe7dXV1XHHFFTQ0NOCc495772233TXXXMPXv/51tm3bBkBSUhJPPfUUt912GzU1NbS0tHD77bczdmz39n7Nud49fqykpMSd6MKr3/zral5cu4dl353dTVWJiMD69esZM2aM12XEhfb2pZktc86VdPRZHboE8tIDVB1oorEl6HUpIiLSzRR0vD/EoKK20eNKRESkuynogNzDY+l05aWIdLPefnooFpzoPlTQET50CRo0LiLdKxAIUFVVpbA7AYfXowsEAse9DV11CeRpdhQR6QEFBQWUlZVRWVnpdSm92uEVxo+Xgg5IDySQkujXoUsR6VaJiYnHvSq2dB8duiS8qF9eRkBBJyIShxR0EXmaHUVEJC4p6CLyNDuKiEhcilrQmVnAzJaY2UozW2tmP2inzXVmVmlmKyK3G6NVX256gIq6BkJaaVxEJK5E82KURuAc51y9mSUC/zKzF5xzb7dpN885d0s7n+9R+RkBmoOOqgNN5KQlR/vrRUSkh0StR+fC6iNPEyO3mOk+5UbG0u3VBSkiInElqufozMxvZiuACuAl59zidppdaWarzOwpMys8ynZuMrNSMyvtrvEpGksnIhKfohp0zrmgc24SUABMN7NxbZr8HShyzk0AFgCPHGU7DzjnSpxzJTk5Od1S2+HZUXarRyciElc8uerSObcfeBW4sM3rVc65wzMr/w6YGq2actKS8fuMverRiYjElWhedZljZgMij1OA84ANbdrkt3p6ObA+WvX5fUZOarIGjYuIxJloXnWZDzxiZn7CAfukc+5ZM7sbKHXOzQduM7PLgRagGrguivWRq7F0IiJxJ2pB55xbBUxu5/W7Wj3+JvDNaNXUVn56gM2V9R03FBGRXkMzo7SSlxHQOToRkTijoGslNz1AXWML9Y0tXpciIiLdREHXSr7G0omIxB0FXSuaHUVEJP4o6FrR7CgiIvFHQdfK4dlRNJZORCR+KOhaSUnyk5GSqB6diEgcUdC1kZceUI9ORCSOKOja0ErjIiLxRUHXhnp0IiLxRUHXRm5GgH31jTQHQ16XIiIi3UBB10Z+RgDnoKKusePGIiIS8xR0bRwZYqDzdCIicUFB14ZmRxERiS8KujYOz46yWz06EZG4oKBrY2C/RJISfOrRiYjECQVdG2YWHmKgHp2ISFxQ0LVDQSciEj8UdO3Iy9CgcRGReKGga8fhoHPOeV2KiIicIAVdO3LTAzS1hNh/sNnrUkRE5AQp6NqRryEGIiJxI2pBZ2YBM1tiZivNbK2Z/aCdNslmNs/MNpvZYjMrilZ9rWnQuIhI/Ihmj64ROMc5NxGYBFxoZjPatLkBeM85dxJwL/CzKNZ3xOFB47ogRUSk94ta0Lmw+sjTxMit7dUeVwCPRB4/BZxrZhalEo8YlJaMmQ5diojEg6ieozMzv5mtACqAl5xzi9s0GQLsBHDOtQA1QFY727nJzErNrLSysrLb60z0+8hOTWavgk5EpNeLatA554LOuUlAATDdzMa1adJe7+1D1/g75x5wzpU450pycnJ6olQtwCoiEic8uerSObcfeBW4sM1bZUAhgJklABlAdVSLi8jL0OwoIiLxIJpXXeaY2YDI4xTgPGBDm2bzgc9GHl8FvOI8GrWtHp2ISHxIiOJ35QOPmJmfcMA+6Zx71szuBkqdc/OBPwB/NLPNhHty10axvg/IywhQc6iZQ01BUpL8XpUhIiInKGpB55xbBUxu5/W7Wj1uAOZGq6ZjObLSeG0Dw7P7e1yNiIgcL82MchRHxtLpPJ2ISK+moDsKzY4iIhIfFHRHkaf5LkVE4oKC7ihSkxNIS05Qj05EpJdT0B1DrsbSiYj0egq6Y8hLD7BbPToRkV5NQXcMeRkBzXcpItLLKeiOIS89QGV9I8GQJ5OziIhIN1DQHUNuRoBgyLGvvtHrUkRE5Dgp6I4hP11DDEREejsF3TFodhQRkd5PQXcMh4NuZ/VBjysREZHjpaA7hqz+SYzOS+Pp5WV4tFqQiIicIAXdMZgZ188azoY9dSzaUuV1OSIichwUdB24fNJgMvsn8eCb27wuRUREjoOCrgOBRD+fOnUoL2+oYNu+A16XIyIiXaSg64RPzRhGgs945K3tXpciIiJdpKDrhEHpAS6bMJgnS3dSc6jZ63JERKQLFHSd9LlZwznYFOQvpTu9LkVERLpAQddJ4wsymF6UyUNvbqclGPK6HBER6SQFXRdcf3oRu/YfYsH6vV6XIiIinRS1oDOzQjNbaGbrzWytmX2lnTYfMbMaM1sRud0Vrfo6Y3ZxHgUDU3jwX9u9LkVERDopmj26FuAO59wYYAbwZTMrbqfdG865SZHb3VGsr0N+n3HdaUUs2V7N6rIar8sREZFOiFrQOed2O+eWRx7XAeuBIdH6/u5y9bRC+if5eUgDyEVEegVPztGZWREwGVjcztszzWylmb1gZmOjWlgnpAcSmVtSyN9XlVNRq1UNRERiXdSDzsxSgaeB251ztW3eXg4Mc85NBP4X+NtRtnGTmZWaWWllZWXPFtyO604roiXk+NPbO6L+3SIi0jVRDTozSyQcco855/7a9n3nXK1zrj7y+Hkg0cyy22n3gHOuxDlXkpOT0+N1t1WU3Z9zR+fyp8Xv0tAcjPr3i4hI50XzqksD/gCsd87dc5Q2eZF2mNn0SH0xuWzA9acXUX2gifkryr0uRUREjiEhit81C/g0sNrMVkRe+xYwFMA5dz9wFfBFM2sBDgHXuhhdCG7miCxG56Xxuze2ctXUAnw+87okERFpR9SCzjn3L+CYaeCcuw+4LzoVnRgz4+azRnL7vBUsWL+X88fmeV2SiIi0QzOjnIBLJ+RTmJnCr17dohXIRURilILuBCT4fXzhzJGs3LlfK5CLiMQoBd0JumpqATlpyfz61S1elyIiIu1Q0J2gQKKfG04fzr8272NV2X6vyxERkTYUdN3gk6cOJT2QwK8XqlcnIhJrFHTdIC2QyGdmFvGPtXvYXFHndTkiItKKgq6bfG5WEYFEH795davXpYiISCsKum6SlZrMtdOG8syKXezaf8jrckREJEJB140+f+YIAH73unp1IiKxQkHXjYYMSOGjk4fwxNJ32Vff6HU5IiKCgq7b3XzWSBpbQlqYVUQkRijoutlJg1K5oDiPRxftoK6h2etyRET6PAVdD/jS2SOpa2jhT2+/63UpIiJ9XjSX6ekzJhQM4IxR2fxq4WZqG5r5xPShFGb287osEZE+ST26HvKjj45j5sgsfvvaFs78r4V87qElvLx+L8GQVjkQEYkm6+3Ly5SUlLjS0lKvyziq8v2HeGLpTp5Y8i4VdY0MGZDCx6cXcvW0QgalBbwuT0Sk1zKzZc65kg7bKeiiozkYYsG6vfxp8Q7e3FxFUoKPBz87jdNHZXtdmohIr9TZoNOhyyhJ9Pu4aHw+j904g1fuOIuirH585Yl32Fvb4HVpIiJxTUHngRE5qfz6k1M41Bzk1j+/Q0sw5HVJIiJxS0HnkZMGpfHTOeNZsr2a/35pk9fliIjELQWdh66YNIRPnDqU37y6hVc27PW6HBGRuKSg89hdlxZTnJ/OV+etpOy9g16XIyISd6IWdGZWaGYLzWy9ma01s6+008bM7JdmttnMVpnZlGjV55VAop9ff3IKoZDjlj+/Q1OLzteJiHSn4wo6M0s6jo+1AHc458YAM4Avm1lxmzYXAaMit5uA3xxPfb1NUXZ//vOqCazYuZ+fvrDe63JEROLK8fbofnD4gZnN6swHnHO7nXPLI4/rgPXAkDbNrgAedWFvAwPMLP84a+xVLhqfz+dmFfHQm9t5YfVur8sREYkbxxt0/zSzm8zsk4R7YV1iZkXAZGBxm7eGADtbPS/jw2FI5LtLzay0srKyq18fs7550RgmFg7g/z21iu37DnhdjohIXOhy0JnZH4DLgWuBsc6573Tx86nA08Dtzrnatm+385EPTd3inHvAOVfinCvJycnpytfHtKQEH7/6xGR8PuOWx5fT2BL0uiQRkV6vy0HnnLsB+DbwXaDCzH7b2c+aWSLhkHvMOffXdpqUAYWtnhcA5V2tsTcrGNiPn8+dyJpdtfz0+Q1elyMi0usd16FL59wh59ybzrlfOOe+0JnPmJkBfwDWO+fuOUqz+cBnIldfzgBqnHN97oTV7OJcrp81nIff2s4/1vS5H19EpFt1aT06MysExgLjgPGED112OKFmxCzg08BqM1sRee1bwFAA59z9wPPAxcBm4CDwua7UF0/uvGg0y3ZU8/WnVjF2cIbWsxMROU4drl5gZl8APgsUA8nAc8AaYDWw2jnn6fxVvWX1guOxs/ogF//yDUbkpPKXL8wkKUHj+0VEDuvO1Qu+CXwVmAo8CwSAB51zT3sdcvGuMLMf/3nlBFbu3M9//VPn60REjkdngu5S59xi59wW59xc4D7g72b2VTNTF6OHXTQ+n8/MHMbv3tjGy+s1H6aISFd1GFTOuTVtnv8DmA5kAm/2UF3SyrcuHkNxfjp3/GUl5fsPeV2OiEivcrxXXTY6575L+Nyd9LBAop9ffXIKzS0hbntc69eJiHTFCR161Dm66Bme3Z+fzBlP6Y73+MbTq6nQyuQiIp3SpeEF4q0rJg1h3e5aHnh9K/NX7uKKSUO48YzhjM5L97o0EZGY1eHwglgXz8MLjmZH1QEe/Nc2niwt41BzkDNGZfP5M0ZwxqhswuPyRUTiX2eHFyjoerH9B5t4bPG7PPzWdirrGhmdl8ZNZ47go5OG4PMp8EQkvino+pDGliDzV5Tz+ze2sXFvHSXDBvLjj43nlLw0r0sTEekx3TlgXGJccoKfuSWF/OP2M/j53Ilsqaznkl++wX+8sIFDTVoBQUT6NgVdHDEzrppawCt3fIQ5U4Zw/2tbmH3vayzcUOF1aSIinlHQxaGB/ZP4z6smMu+mGQQS/Xzu4aV86bFl7NWQBBHpg3SOLs41tYR44PUt/O8rm/GZkZuejN9nJPp9JPiNBJ+PxMj9uWMGcd1pRST49f8fEYl9uhhFPmBH1QF+/8Y2ahuaaQk6moMhWkKR+6Cj5lAz63bXMnZwOv8xZwLjCzK8LllE5JgUdNIlzjleWLOH789fy776Rq47bTh3nH8y/ZM1p4CIxCZddSldYmZcPD6fBXecxSdOHcpDb21j9j2vsWCdVkwQkd5NQScfkB5I5EcfHc9TN59GWiCRGx8t5Yt/WsbuGq2aICK9kw5dylE1tYT43Rtb+Z+X/00w5Jh1UjaXTcjngnF5pAcSvS5PRPo4naOTbvNu1UGeWPou81eWU/beIZISfJx9Sg6XTRzMuaNzSUnye12iiPRBCjrpds45Vuzcz/yV5Ty3ajcVdY30S/Jz8fh8vndZMWnq5YlIFHU26HRJnXSamTF56EAmDx3Idy4pZvG2Kv6+spy/lJaxaW8dj14/nQH9krwuU0TkA3QxihwXv884bWQ2P50zgd9+eiob9tRx7QNvU1nX6HVpIiIfELWgM7MHzazCzNYc5f2PmFmNma2I3O6KVm1yYs4dk8uDn53GjqqDXPPAIl2hKSIxJZo9uoeBCzto84ZzblLkdncUapJucvqobB69YToVtY3MvX8R71Yd9LokEREgikHnnHsdqI7W90n0TSvK5M+fP5X6xhau/u0iNlfUe12SiEjMnaObaWYrzewFMxt7tEZmdpOZlZpZaWVlZTTrkw5MKBjAEzfNoCUU4toHFrF+d63XJYlIHxdLQbccGOacmwj8L/C3ozV0zj3gnCtxzpXk5ORErUDpnNF56cz7wkwSfD6ufeBtXtI0YiLioZgJOudcrXOuPvL4eSDRzLI9LkuO08icVP5y80xy05P5/KOlXP/wUnZUHfC6LBHpg2Im6Mwsz8ws8ng64dqqvK1KTkRhZj+eu+0Mvn3xGBZvrWL2va9zz4sbOdQU9Lo0EelDojZg3MweBz4CZJtZGfA9IBHAOXc/cBXwRTNrAQ4B17rePm2LkOj38fkzR3D5pMH8+Ln1/PKVzTy9fBd3XVbM+cW5RP5vIyLSYzQFmETVoi1VfG/+Gjbtreesk3OYM2UIfp+R4DN8ZiT4I/c+H4WZKQzL6u91ySISozTXpcSs5mCIRxft4N6XNlHf2HLUdn6fcePpw7n9vJM1cbSIfIiCTmJeXUMzFXWNhEKOlpAjGLkdfvx/75Tx+JKdDM3sx08+Np7TR+naJBF5nyZ1lpiXFkg85ooH04dncsWkIXzrr6v51B8WM2fKEL57STED+2viaBHpvJi56lKkPTNGZPH8V87glrNPYv6Kcs695zWeWbGL3n4kQkSiR4cupdfYsKeWO59ezYqd+zn9pGwunZDPtOGZjMjur6s3RfogHbqUuDM6L52nv3gaf1y0nfsWbuZfm/cBkJ2aRMmwTEqKBjJ9eCbF+ekk+HWwQkTC1KOTXsk5x9Z9B1i6rZol26tZur2andXh5YFSkxP4xoWn8OmZRd4WKSI9Sj06iWtmxsicVEbmpHLt9KEA7KlpYMn2av5SupPvPrOWyvomvnreKB3WFOnjdHxH4kZeRoDLJw7moeumMXdqAb98+d98529rCIZ691ELETkx6tFJ3Enw+/jPqyaQlZrM/a9tYf/BZu65ZiLJCRp0LtIXKegkLpkZd140mqz+Sfz4+fXsP9TEbz9dQmqy/smL9DU6dClx7fNnjuDncyfy9tZqPvG7t6mqb/S6JBGJMgWdxL2rphbw209NZeOeOubev4h/763zuiQRiSINL5A+Y+n2aq5/eCl1DS0UDEzhtJFZzDopm5kjsxiUFvC6PBHpIk3qLNKO3TWHeHHtXt7aso9FW6qobQivnjBqUCqzTspmyrCBDBkQIDc9wKC0AEkJOughEqsUdCIdCIYca8treGtLFW9u3sfS7dU0NIc+0CY7NZm8jGTy0gMUDOzH7OJcZozIwu/T2DwRrynoRLqosSXI1soD7KltYG9NQ/i+toHdNQ3sqWng3eqDHGwKkp8R4KOThzBn8hBG5aZ5XbZIn6WgE+lmDc1BFqzfy1+X7+K1TZUEQ44JBRnMmTyEyyYOJis12esSRfoUBZ1ID6qsa2T+ynL+uryMteW1+H3GwH6J+H1Ggs8XuTd8kfvs1GSKB6dTnJ/O2MHpDM/ur4mnRU6Qgk4kSjbuqeO5VeVUHWgi5BwtwQ+ulN4cDLGntoENe+poagmfAwwk+jglLxx8Z52cw4Xj8jz+KUR6H03qLBIlp+SlcUreKR22aw6G2Fp5gLXlNawtr2VdeS3PrSrn8SXv8ucbT+W0k7KjUK1I36MenYiHGpqDnHfPa/RPSuC5207X4UyRLuhsjy5qv1Vm9qCZVZjZmqO8b2b2SzPbbGarzGxKtGoT8Uog0c93Lilm4946Hlv8rtfliMSlaP738WHgwmO8fxEwKnK7CfhNFGoS8dwFY3OZdVIW//3iRqoPNHldjkjciVrQOedeB6qP0eQK4FEX9jYwwMzyo1OdiHfMjO9dNpYDTUH++8WNXpcjEndi6YTAEGBnq+dlkdc+xMxuMrNSMyutrKyMSnEiPenk3DQ+PWMYjy95l7XlNV6XIxJXYino2ptTqd0rZZxzDzjnSpxzJTk5OT1clkh0fPW8k8lISeQH89fR2y8SE4klsRR0ZUBhq+cFQLlHtYhEXUa/RL5+wWiWbK/m2VW7vS5HJG7EUtDNBz4TufpyBlDjnNNvu/Qp10wrZOzgdH7y/HoONrV4XY5IXIjm8ILHgUXAKWZWZmY3mNnNZnZzpMnzwFZgM/A74EvRqk0kVvh9xvcvH8vumgbuf3WL1+WIxIWozYzinPt4B+874MtRKkckZk0ryuTyiYO5//WtzC0ppDCzn9clifRqsXToUkQivnnxaPxm/Pi59V6XItLraa5LkRiUn5HCl88eyc9f3MSFv3id4sHpjB2cwdjB6RQPTic9kOh1iSK9hoJOJEZ9/swRACzd/h6vb9rHX5fvOvLe0Mx+jB2czuShA5g+PIuxg9NJ1DyZIu1S0InEqOQEP7ecM+rI84raBtaW1x5Z/WBNeQ0vrNkDQL8kP1OHDWR6USbTh2cysXAAgUS/V6WLxBQFnUgvMSg9wKD0AGePHnTktYraBpZsr2bptmoWb6vmngWbcA6S/D6mDhvI+WNzmV2cS8FAXdAifZeW6RGJI/sPNlG6/T2WbK/m1Y0VbNpbD8DYwemcX5zHBeNyOSU3DbP2JiIS6V20wriIsG3fAV5cu4cX1+1l+bvv4Vz4/N6ZJ2eT1T+Z1OQEUgMJR+7TkhNIT0lkeHZ/nfOTmKegE5EPqKhrYMG6Cl5ct4dl29+jrvHoM6/0T/Jz6ogsThuZxWkjsxmdl4bPp16gxBYFnYgcUyjkONDUQn1jCwcaW6hrCD+uPtDE0u3VvLW5iq37DgAwsF8iMyOhVzAwBeDI4c/D8WcG/ZISKByYQk5asg6PSo/rbNDpYhSRPsrnM9ICiaS1MybviknhFbJ21xzirc1VvLWlire27OP51Xs6te3kBB8FA1MozOxH4cB+FGamMHhACoEEP4kJPhL9RpLfR4I//DiQ6Gd4Vn/1GqVHKOhE5KjyM1K4cmoBV04twDnHjqqDVB9s4v0DQeEHh5/XNbSw872D7Kw+yM7qQ+x87yDLd7xHbUPHE1SPG5LOnReO4fRR2T3zw0ifpaATkU4xM4qy+1NE/y5/tuZQM3trG2hqCdEUDNHcEqI56GgKBmlqcVTWNXD/a1v51B8Wc8aobO68aDRjB2f0wE8hfZHO0YlITGhsCfLHRTu4b+Fmag4189FJQ7jj/JM1BlCOShejiEivVHOomd+8uoWH3tyGc/CZmcP44kdGkpWa7HVpEmMUdCLSq5XvP8S9L23iqeVlJPiMc0YP4sopBZw9epDG+AmgoBOROLG5oo4nluzkbyt2sa++iaz+SVw+aTBXTilg7OB0DWPowxR0IhJXmoMhXt9UydPLy1iwroKmYIjReWl8bPIQLhibR1F21y+Skd5NQScicWv/wSb+vrKcp5bvYuXO/QCcnJvK+cV5nD82l/FDMtTT6wMUdCLSJ+ysPshL6/by4ro9LNlWTchBXnqA2cW5fOSUHDJSEvH5jASf4Y/cwo995KSF5/uU3klBJyJ9znsHmnhlQ3g+z9c37eNQc7DDz2T2T2JoZr8P3Aoz+1Gcn05GP63kHssUdCLSpzU0B1m9q4bG5hAtoRDBkHv/5hzNwRB7axvZURWeyeXd6oPs2n+IYCj8NzHJ7+O84kFcNbWAM0flkKArPWOO5roUkT4tkOhnWlFmlz7TEgyxu6aB7VUHWLihkr+t2MXzq/eQk5bMxyYP4aqpBZycm9ZDFUtPiWqPzswuBP4H8AO/d879RztDucYAAAqASURBVJv3rwP+C9gVeek+59zvj7VN9ehEpKc0tYRYuLGCp5aVsXBDBS0hx4SCDGaPySU50YdF1m5ofd2L2fvnA4/c+8PnBBN8xoic/pw8SMsedYeYO3RpZn5gEzAbKAOWAh93zq1r1eY6oMQ5d0tnt6ugE5Fo2FffyDMrynlqWRnrd9ee0LYG9ktk+vBMTh2exakjMhmTl67gOw6xeOhyOrDZObcVwMyeAK4A1h3zUyIiMSA7NZkbTh/ODacP52BTC5FTeTjnONxdcA5wEHSOlmCIlsg5wfB9iMaWEOvKa3l7azWLt1Xxz7V7AUgPJBwJvmnDMxk7OF2zv3SjaAbdEGBnq+dlwKnttLvSzM4k3Pv7qnNuZ9sGZnYTcBPA0KFDe6BUEZGj65d0/H86xw7OYG5JIQC79h9i8dYqFkeCb8H6isj2/UwZOpBpRZlMH57J5KEDCCT6u6X2viiahy7nAhc4526MPP80MN05d2urNllAvXOu0cxuBq52zp1zrO3q0KWIxIuK2gaWbK9m6bZqFm+rZuPeOpyDRL8xZehAbjpzBOeMHqTB8BGxeOiyDChs9bwAKG/dwDlX1erp74CfRaEuEZGYMCg9wKUTBnPphMFAeCWHZTuqWbLtPZ5bXc4Nj5Qybkg6t50zitnFuQq8Topmjy6B8OHIcwlfVbkU+IRzbm2rNvnOud2Rxx8DvuGcm3Gs7apHJyJ9QXMwxN/e2cV9Czezo+ogxfnp3HbuKM4vzu2zF7LEXI/OOddiZrcA/yQ8vOBB59xaM7sbKHXOzQduM7PLgRagGrguWvWJiMSyRL+PuSWFfGzyEJ5ZUc59Czdz85+WMTovjS+ffRKjclNJ8BkJPh8J/vfvE30+0lMS+nTvTzOjiIj0Qi3BEH9fVc7/vrKZrZUHjtm2YGAKcyYP4WNTChgeR6s8xNw4up6ioBORviwYcry9tYraQ800R4YxNAcdLcHw44bmEK//u5J/bd6HczBl6ADmTCngsgmDe/1cngo6ERE5Yk9NA8+s2MXTy8vYtLf+yFyel08czOmjcnrlKg4KOhER+RDnHGvLa3l6eRnzV5RTdaCJJL+P6cMzOXv0IM4ZPajXHN5U0ImIyDE1B0OUbn+PhRsreGVDBZsr6gEoyurH2aMHMbs4l5kjsmL2QhYFnYiIdMnO6oO8sqGChRsreGtLFU0tIc4bM4i7rxjH4AEpXpf3IQo6ERE5boeagjy2eAf//eImfAbfuGg0nzp1WEyN2ets0GnWUBER+ZCUJD83njGCF796JlOGDeSuZ9Zy1f1vsWlvndeldZmCTkREjqowsx+PXj+de66eyNZ9B7jkl29w70ubaGwJel1ap/W+60lFRCSqzIw5Uwo46+QcfvjsOv7n5X/z91XlTBk6EJ+BYfh8AIbPwGfG+WNzOWNUjtelAzpHJyIiXbRwYwX/9Y+N7D/YRMiBw4XvXXj4QkNzkIaWEPd9fDIXjc/vsTpibq5LERGJD2efMoizTxl01PfrG1v47INLuPXxd7jP4MJxPRd2naFzdCIi0q1SkxN45PrpTCwcwC1/fod/rNnjaT0KOhER6XapyQk8/LlpTCjI4JY/L/c07BR0IiLSI9ICiTxy/XTPw05BJyIiPeZw2I2PhN0/10Y/7BR0IiLSow6H3bghGXz5seW8GOWwU9CJiEiPSw8k8ugN4bD72l9WUnOoOWrfreEFIiISFYfDbnNFPRkp0Vv0VT06ERGJmvRAIlOGDozqdyroREQkrinoREQkrkU16MzsQjPbaGabzezOdt5PNrN5kfcXm1lRNOsTEZH4E7WgMzM/8CvgIqAY+LiZFbdpdgPwnnPuJOBe4GfRqk9EROJTNHt004HNzrmtzrkm4AngijZtrgAeiTx+CjjXzGJnOVsREel1ohl0Q4CdrZ6XRV5rt41zrgWoAbLabsjMbjKzUjMrrays7KFyRUQkHkQz6NrrmbVdDK8zbXDOPeCcK3HOleTkxMbCfiIiEpuiGXRlQGGr5wVA+dHamFkCkAFUR6U6ERGJS9EMuqXAKDMbbmZJwLXA/DZt5gOfjTy+CnjF9fYl0EVExFMWzRwxs4uBXwB+4EHn3I/N7G6g1Dk338wCwB+ByYR7ctc657Z2sM1KYMcJlpYN7DvBbcQ77aOOaR8dm/ZPx7SPOtZ6Hw1zznV4/iqqQRerzKzUOVfidR2xTPuoY9pHx6b90zHto44dzz7SzCgiIhLXFHQiIhLXFHRhD3hdQC+gfdQx7aNj0/7pmPZRx7q8j3SOTkRE4pp6dCIiEtcUdCIiEtf6fNB1tHRQX2RmD5pZhZmtafVappm9ZGb/jtxHd4ngGGJmhWa20MzWm9laM/tK5HXtowgzC5jZEjNbGdlHP4i8PjyyBNe/I0tyJXldq5fMzG9m75jZs5Hn2j+tmNl2M1ttZivMrDTyWpd/z/p00HVy6aC+6GHgwjav3Qm87JwbBbwced5XtQB3OOfGADOAL0f+3Wgfva8ROMc5NxGYBFxoZjMIL711b2QfvUd4aa6+7CvA+lbPtX8+7Gzn3KRWY+e6/HvWp4OOzi0d1Oc4517nw3OMtl5C6RHgo1EtKoY453Y755ZHHtcR/kM1BO2jI1xYfeRpYuTmgHMIL8EFfXwfmVkBcAnw+8hzQ/unM7r8e9bXg64zSwdJWK5zbjeE/9ADgzyuJyaYWRHhKesWo330AZHDciuACuAlYAuwP7IEF+j37RfA/wNCkedZaP+05YAXzWyZmd0Uea3Lv2cJPVhgb9CpZYFE2mNmqcDTwO3OuVqtEfxBzrkgMMnMBgD/B4xpr1l0q4oNZnYpUOGcW2ZmHzn8cjtN++T+aWWWc67czAYBL5nZhuPZSF/v0XVm6SAJ22tm+QCR+wqP6/GUmSUSDrnHnHN/jbysfdQO59x+4FXC5zMHRJbggr79+zYLuNzMthM+ZXIO4R6e9k8rzrnyyH0F4f8sTec4fs/6etB1ZukgCWu9hNJngWc8rMVTkXMpfwDWO+fuafWW9lGEmeVEenKYWQpwHuFzmQsJL8EFfXgfOee+6ZwrcM4VEf6784pz7pNo/xxhZv3NLO3wY+B8YA3H8XvW52dGaW/pII9L8pyZPQ58hPByGHuB7wF/A54EhgLvAnOdc31yUVwzOx14A1jN++dXvkX4PJ32EWBmEwhfKOAn/B/qJ51zd5vZCMI9mEzgHeBTzrlG7yr1XuTQ5decc5dq/7wvsi/+L/I0AfhzZGm3LLr4e9bng05EROJbXz90KSIicU5BJyIicU1BJyIicU1BJyIicU1BJyIicU1BJyIicU1BJyIicU1BJ9ILmdl5ZvZHr+sQ6Q0UdCK900TCM2eISAcUdCK900TgHTNLNrOHzewnpuUTRNrV15fpEemtJhKetf2fwO+dc3/yuB6RmKW5LkV6mcgSQfuAHcAXnHOLPC5JJKbp0KVI71NMeImpFiDocS0iMU9BJ9L7TATeIryO2UNmlutxPSIxTUEn0vtMBNY45zYB3wCejBzOFJF26BydiIjENfXoREQkrinoREQkrinoREQkrinoREQkrinoREQkrinoREQkrinoREQkrv1/BwULINhBFeIAAAAASUVORK5CYII=\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 \"*Multiple Testing and Variable Selection along Least Angle Regression's path*\" (Azaïs and De Castro, 2019), 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 \"*Fast algorithms for component-by-component construction of rank-1 lattice rules in shift-invariant reproducing kernel Hilbert spaces*\" by D. Nuyens, R. Cools iun Mathematics of Computation, 2006. 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.5104805749849506\n"
     ]
    }
   ],
   "source": [
    "from 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 5000 Monte Carlo simulations and 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 tyour 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.4"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": []
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}