Parameter fitting (Least square and Least absolute deviation) with different error distributions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We numerically explore the relationship between error distributions, objective functions (least squares versus least absolute deviation) and computational speed.\n",
    "\n",
    "The following simulation creates a data set according to\n",
    "$$y_i = \\beta_1 + \\beta_2 x_i + \\beta_3 x_i^2 + \\epsilon_i.$$\n",
    "\n",
    "We generate $n$ observations from this model for $x$ on an equally spaced grid $[0,1]$. For $Yt$, the errors ($\\epsilon_i$) have a t-distribution with 2 degrees of freedom. For $Yn$, the errors have a standard normal distribution. We plot each of these data sets and the true regression function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "set.seed(1234)\n",
    "n <- 500\n",
    "x <- seq(from=0,to=1,length.out=n)\n",
    "X <- cbind(1,x,x^2)\n",
    "beta <- matrix(c(1,2,4),nrow=3)\n",
    "Yt <- X%*%beta + rt(n,df=2)\n",
    "Yn <- X%*%beta + rnorm(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4QAAAGkCAMAAABeoxsWAAAAM1BMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD/AAD///89ODILAAAACXBIWXMAABJ0\nAAASdAHeZh94AAAgAElEQVR4nO1diWKrOAxUer9um/r/v3aTgKSRbYghJtiJZvflAF94PDoM\nbSk4HI5dQXsPwOF4drgIHY6d4SJ0OHaGi9Dh2BkPLsL3vQfg2B5Ena/izoc/j/8OD315jgEu\nwpbRPTuOEnRPc+fDn0f37DhK0D3NnQ9/DjQAjhw/DnT4+OWzvy/0AR9C+H4/fXn/DrbA8fP1\n9Pnt373H/+g4c/N9mtr3n/FAfvrPxb5e6OW/EL4O9PrfcPrf2+n4y8BlJML+WH4mEf4ehiP/\njWdfiF7hQ3gda7wFLMC1LkUc9XCa0o9hZgcV5qf//GE48/uh7HHZyzcrwg5ZfiYR8jwf9Ow/\n+PDGxA788PH3y8njibuvfa7jUSHTTZc97Inp12IHPft1EssxhI/hmxVhhyw/sAhjdgbmjsN0\nXwzsMcCH79OHr9P5z9P7Nxw/2+DT2/FkMu9/BY+Ms1BOE/0+sDQ1/eOHr7PL+rm8nY69DJyM\nBBuae2T5iUR4soEjH6P5/OZClw/vbAQ/BtPMx8+WlTMIR0WME3wcWJqa/vOH/8ybbSI+1iPL\nTyRCDX8Ow7cjHz7ie/gVai/fP8eQyXVYGcyOne54+rmYfTsV+/dxye5CRHOPLD+lCA1zdi3o\nJ/n+waT+3mnYT4JEhOZ4zI99+/eSo3Is0B3LTyTCQyaSCYYesZEHU/n4b9hR8+3Rqpj2hGb6\ncyL8d3p7ef/6yYiwR5afSIRvnB/Yc/zhLc0WoPL3u/3uuBlWhJPTnxPhy0hlRoQ9svzQC+s8\no0f5djKfh/8ub68hQ09m3+xy/EWyicNdB//wsCKcmv6sCMcvOU/YI8sPLcLzlteH/XqB2WaT\nD3IDeLhvBbS9/l4y94/gqIgoTZuY/qwIXy9kfB8yIuyR5YcW4bsN8b/H2b9Mc0qP8PNuj3PK\n7ilhXUQinJr+nAj/G4seLlKLQsr+WH5oEZ4zAPyJwuPHKep4kxtHwX44pwQH81ShHD6T48/L\nVEYswonpz4kw/JwoObz//F7uBloRdsjyY4vQ4egALkKHY2dsLkJyVMLWTDmBW2NyildS89/n\n8Dj628d/Vzhc2YEjQu2JLGXQCayFyiI8voC+5/eTnMNJzNjGXOmqfZcz6ATWQmURftDh3/CT\nmL/fh/k7K87hBMZbXOXlq/ZezqATWAuVRXigH/n8M/+MgXM4AYLX8vK1UM6gE1gLlUUYPSS7\nqucnB0XvxRUqdV/MoBNYC+4JW8POInRPeH/Uzwm/h5+78pxwJXYWYTmDTmAt1L5F8Qp7ay/H\nuZLO4QT2zQnLGXQCa6H+fcKPy12mw9un3ydch313R8sZdAJroboIb+/Zsed9wvb7fTzcV4RF\nz+o4luGeE+kEbgH3hC2jaKW7J+wdLsJ2UZgcugh7h4uwXRRuk7oIe0f1J2aKswbn8Ar4B8av\nTVTtJ2ZKGXQCa6GyCL9chDFW72DQZEBqm6w7keUMPguB26N2OPpzKP2FOM/B4dKbflhVnOF8\nk5UnspjB5yDwHqieE/6U/mq45+Bw6eMv6ObG32AUV4+brD2RpQw+B4H3QP2NmS94AnhVz48E\nfRC0MCpFNyefTc3k2dLqE1nI4FMQeBf47uimME6t/KYfcS1xhpkm4wP3xlMQeBe4CDeFVcz1\nS44UlqvmInw4uAi3hZVR8T0/cYUhdaBb54SleA4C7wEX4bawMloqwuwNDiJ799BF2DtchFvj\nrJdiERbFrRcJqgpdhL3DRXgXFN+pyASgiTOMGnMR9g4X4V2w4J49xRKMa8Zu1UXYO1yEd8La\np9dSH+oi7BvpSnhEET7Sz5tmlOYi7Bm5mOjxRHjD05oNIqc0zwk7xtyt34nCG2IzEW7a+kSf\nm4k+K0JrZlyEHSEbvjycCO8bpF3Ut6nvzdoUI3oXYUdwEVbva1Dfpr6X+6j/R+xuhYtwBVyE\n2/S1dY+swCkZugh7gueEW/Qk4tiwy9lrchH2BN8drdwTdBi2nLF5X+si7At+n7BqP+Z9Y0eI\nj6Dmh3FvuAhr4RFFeDeM6tvc99Kse3cR9g4X4Q3Q3z+xse8df8DeRfiYcBHehDtFvqz2bGcu\nwvZw5ZdFPvzGTO/I8ce/J8pF2AU0eciSmaQWLsK2kE/+zBZQ/Ft5XYStgXfq5sh0EbaLCY+n\nxPEvt/BnR5uFmMzstnnmfpOLsCXkf9MonggE/4/n7jlCHNNO/TaN3C8zySeALsJGMfE7t4dT\nco6sCl2EzSD/a71chD2BjejUbfkQxt+zRuQibBF84xjfJuJRzwl3wrUbGqysNJfHv08Rh6wu\nwlag+2Zj3q5HsZTvju6H2adeJN/LFIr/PgWJFodDG433Gp6OwKtQ/kSHWcb9PmEG97nlnreL\nQ/9BndvU870sQgrRnw50EbaCKJKh0pXlIrzXj13MRI/AWWYo5k49Bb9P2C748UL4Qz6Wqvwy\ncxHOeKiy6oX6nRYhnsn/2nuIPxek9RujGQLbgTwlQ2NgY6zqpLV3Ed62v1HuRstEOFlRRFie\n1m+MVghsCfwcxZi5W5M5ae1dhDeKsLzuVQ4m9mziisVp/cZohcCmwOoLGLxocBpS+oKL8EYR\nLqk87TRz8hSu5n+YsLDvLdAKgU1hDEUhgyCTxkuoaipNtrblUO/SQSluyQmXKXgqfUxUhpvb\nmmZcGcO90QyBbYG5Y874mz51OPEXfnJtbTnQu3RQilt2R2+LZXEMRoKS2Gvb0z24CNuC/gwT\nRc5wKqxxEYZpD1VSFV5rgcR6cgQjXnHR3trGaIjAlkCsuSDKk5v3Eoq6CGtig5uMxEGLNJyq\nMSq/Cx6DwOrA9IH0uQr7kIVl0UV4K6o/bsMbbEaE048iugjbAns6Yt2J+IgdY4gyfBdhOzDP\nWbAMA3A5Four3XeUu/fbOob9F0kASQGhqouwSWiwyWk8P6ythnUoGFe880D37rd1wP0JTuTN\njQt4uJtrTDa1+Vi37qAzaLCJLCGXE/Goi7A1oNsD1alljbJ7F2EjwJsd5iaTiDC7sZY7cCc4\ngbPQPRjQoPGCjyTC6jsjuwBEqI6PUhESxRfsItwbEysQRAhaBBcYv6ctVB3mdh1scI9guqsN\nu4lFqN2ZtDBzt9dFuCPgDgQe03eKgOng44iwZmOzHW2sdsgJpUe8wwRPYgRzwS7C3RA9XRgf\ni0JRc88C6/QuwqvXUQ1bqz0hVNN6JlFi1L8/GIiLcDcMJtOujOgbPjLD/10YLPkFJS7C+3Qk\nwQkHNsHqUTdqRkLDiUAXYRNgi4luTUMVsJySR4yb3WcGIblxEdbpaGW+KHKLN62BSwlhBsN7\nkWBwEe4Pkru6IfCL5VJLQD5IowQLGOxEhHfLCWdFuDpfxK2XqHE1sxqKjm6QTEkX4S6A7WtN\n5EV1vIsGESuHOqMVLTGj3Yhw4/0S7Qhel5y73iaMH+4/iAhhT2Ywony3Imrk3nh6EdoNM3wh\nkwsGeOVAJjyUCFdHgsu7CVNqXx0TpyKETtiG6rba38Cg3yfcHxJtyqZZQBGyR9RUHhgkTu6l\nscletr2IHjmcVHsiwlK7kPGE8DpSOv4f/jgbpHwj90Z/BFaEuXnEQoOYVDN62WEDCcIm3FBs\nqpftLuBOHdwPkQjnI+TM9Cc5oTQEEAKner83HojAxTC/T83oMLAvDEFc5VhlzCWs9xxOTXaz\nMR6Jw5yEJrJHZWj4GoA06xAh2CHZj8m06iK8Pzjvs3fmbUAagrntNN5Y4rjGRVg7kTSb0jMZ\noqQH+I3AampV2yQng9lGK17IEjy7CPHW7XBQHGIAjY5nYUuN+KcNobXJbrbFnhyu2lKd1a3Y\nRLB9E/1yJAOazeSEBOVHN8jxadzossuohocW4TUjTcwd7FRHfhHjziGX0O9kenhSES4fQZFu\n1Tpmm7f5ujabGFPVoBAIZjYahYuwOvJkmwgSgk8Tv6RFIJ2HSAhvMz2lCGcCxit1rplH9F8T\n5wcTiRGpKU1Qm+DxCn1mJm7aRVgdMM3KUCRMSSbUnhrZir2VdJ6TD5KQ1BSdHMeG6EuERVUI\nlKXcmck2rg6i17j0eJCGXBDDVqIoM3QR1oaSDcqL7V/s/xKRBolECfUMDjFuanIgm+FhRYiJ\nW8QM/0QgizCIRRSragJU2Y8RvWYCJRdhbYAI5Xu6AMTVKXfm9KUW74miozQ7B7bSxEA2w54c\nlsSWuRrX41Hef853hBtq/EAhf9FghYsOBII08/Goi7A2hD30dtE5oMRukeJpedie1IIO58oY\nfGwRpv7kahV4nWsXczfDjOyWjZRqTgGb1jCuP9YgNw3ZhxmHi7A68Onr4XvWCpvdGUv5+cAf\n7mqHRHrY3eQ4Nsa+HHIqVqzEMt2KwbOa1Uew+asUBtfI6j3/42fUopasQwx46u54ZBFG0wwB\nSlqWS8BX3lH7s5va2du82UavnaiFBjhc5hCLBGs3WqLkgkxyoFvYsmMmoac8JhplgfykhYtw\nU0g4I6nB5GKRHB8ZJ7kzGGJmpRq09NwihNfFdfOSjKIW9HzMhlhGSOnNEZVgUJJNt2UZxcZo\ngMCNkAQww7eEctlli1TGDJL8j4FMCLGgn1mEuTC/tOqkD41coCQMSRovESk/ysQF84EoF453\n1lyE9YELA2Y7kWCw6pOSeG8X0ow4dHURhnCTCOVVreF4JpInejrM5A0hun2W+aFPqavO1EVY\nHYZD+z4lQpGgKpEoQCrBBlaz/xD4wEzDpottsT+HS0RoF77RB6YEmaIqoRBUr3K3UNR5gW6o\nGbOpzpOQxPLhb4H9CayHbHwYrw+ccwhuwDqGYLZEgzJnHsuQOChtOBpWhUubRQMcxvHedMHY\nvZnjsqWZrSgxiQ1t5GaGiHDcEiXVaYDWuQ8MZu1g7o0GCKyGOD6MNBkl4qy+gOoLyiAU4EiU\ncAHgP9NvflgbogEOpzO7uCC84qfcjMYV1YNBb5DeqT0dbwxKKVWfyD3XnYvwjImNssLK0btp\nDhJxdWYY2wwssQatODHyhJDUGtKnFmEpdRmOxtohK0JpFc5p8heAV3kj+KHPADSLW5wKZVyE\nZxSb06nqM0EREyZWlMm1tlV/6BNCGBPWqHajlOLJRViGnKEMGKUYEca+jEY6gsgK04jxCNwY\nhK4gfJW1ECX1LsIzZjRUUFmS+lwDLB2deBQhW1H4TYbmcW0T1AQ4UJTVP74IywOYVIRWgbnw\nUBgybm88IfVAghjvJKGKqji12S7C68nVtdpRfGiiUS4ibJEWkrjEbMco6WiSUYR26T2tCBcF\nMJOGFvZcbFnrEA11vLciabr+JsOY4KSZ3MBdhLeJEEIUPoJTTKZI0HJgNf9iK8r2d9SjHWRi\n/Z9XhEsGMaNYgmk2TZKqTCE+TY0lP6iNhtV0mHJWaEc3xiOJ0E5qJBoJazQ9DKJbgi1RlRel\n62LOlM8PbUPsy+FS2hbHrpruycOFSiUbVGICtQcCouMu84NwES40qdm6UFuiF3mPsngMWfSn\nXdBT6vYLdjRlyl2Ey6te0SMvCOJX2cYGYocXvDFoIlF2lrbfkOXQRTgbqhRUhtfxk4iITSJ3\nwjaVS7MXVL4xhonnaGLpuAiXVpykO8riNMoB60nKrvxCZg1FsSGimMMpa+8iPKMkVJnSQEQp\nxyCitSitCBJojl6Qm1DzSlNk5ce1+EQJvl6I3r5X9nwfrA1gJqWAFlBTdOGNdSglOBBFyQat\nLtTajjNdbzKRBQy2JcLriLSGeUKcvA0aNHMullE3SflnPo0IIVcsds2VRTh0+zpc18e6nqsj\nawBXBjBXpIAUgwsMIVgvKL8+JtOQPOmWc4Sbi7Ccwe5ECK8SWmbXgW68mKrIDckPnAXdP+Wo\nVM3ukpEtOTHf3LnaB30cQ/j9oK9VPVfGdDpcJMHsDlc6enuc1ByajJ4PciCaG5ghOuohS+wG\nIixisDMRRgyFwIIJ6VKQLFC+ByX1chx++QhkD8SbMhLMpMNYktXfIMIDHc+fj/QyW3RVB8vB\nk16mubhyLJQSEQohAaWnoY/86goIg0CixOFP1MdUiLOBCIsY7EKEMMP4Ls5LRZSxtYSLRx8w\nDDGDEPXAkZzlz1O4hQi5k5yVB6zqYPmIMgNbWpuiI5nRGxEGFJL4xpEmuS/I5XkyNFxVMk0X\nOZsd1jI1hXkGdyDwFuDFzIgwqN4CS2o4qQ1whBl0Qw2qwVmCb+mAzFDs0cniS3Hp+Z0pPMwW\nXdXBYpjltLTPjN+bWoMY5UQkc6XADGJ5WCc0DjRjmblsLsjZQIRFDPYgQnkVwsAmBuGHrWYI\nID+VV8BC5k/OK7HS42hs87Fo9J7/WnBiFkRvn1/f9O/08fjRRF4fBxYrKkfV8mGhpANIKqfp\nLF39HWpKtoxOHvzVNXJ9MMuv6QrKGWxfhDJj6NDkAAk7bDQhcMk1df6Ht5XQhqoWkdjJAeUP\nT5ZfBvASRIfjbNFVHawY0uWF4MvCulE1NYKmpAaQ4tBYSzwn0ROGwwsnEbEIMzLPDGbFNV1B\nOYM9iXB802UQsTA6S3CRXAOaIrOhJkqEquJ0JzZG7yPC8PPz9fX2dq58+JjV4P1EiH5raaep\nXRTSIhGmr2T0BrfmMeMDu5neGr46mFXXdAXFDHYkwty7GBuyEgzwFszKsc+I8kuQQhLx5MjD\nERSb0c2n+H4czkQZ16uGYOVmww45iu822wgQiEIwFIT5wPqTMVKu4/yxNRdVCe2LENM8/Zq8\nmzh1XC2WvvObeEHO2jUDgfwRnCl3Yr+EDIXPIMJLb/kFXFIx9ng5P6Vq0hRQOcYH1GxwZG0u\niBcM7PRgZEi7oAcRWqoy7+oG1T2GwSoKTecXecpXEYBDUXvsBpN1l6PwWUQ4HSEsbEXiGGgv\nyitAThfwE4aSMwRtxPLI0Y04xZIhuwhnoFEHDlimV/kwMSqkG6SBqGGLOCdk2xqA2ZD0ND/I\nxSdqoVIHdbRV3t3Yp1IbJIKUsJJw6vUnJURiIXAaaMww9ELabNmQdkAXIhxgszsNP5Snca6t\nCMcP8owoQRg6is8QC/ljsJ/m56pzEa6OMtf3GMQAQrzJYQieDbKdBs6TDIOSXExdlovwZpgY\nA+MQKxHmhd84pTBOELgT2jGkQTMb8NNji7BeU7PdYOgpFlHcIh8Fhzickp82C8EIUZJGIkIR\nUvQpc2Gp53cRziDNyZJXzjEC6CnwJ7Si+MpK5T6UUDMxzyDCsmu8uRcKVoaQLUgIGjAevfwX\nPaWN0Y/GQ3C3AvyfBkzJUCg5/MwiTE1SdAzvwA/foQiKMGSyc/1pMz4lYgtEQCgzb00q9Hfl\nKhafqIWORAgyhDhGiWThIbvwgCEft9GPsB5dhfYzKcHnE2FObDqB9hi/aD2YTGPnxJqqmjg+\nIYr+wCAZyrgBrU7I3tRwpq5u8YlauLUDXKabDhZEBw5MIxB1lWou8Sd2sSF5l/BHKYq0mBsJ\nkJ42e2/cr9/JsCBkVj28yozpNIvNNCFK0Illx5e5J5GDVDfhjB3l1YnqVYQSld/e1NWuwE5K\nd1HESJrCjfz9KenaEL6TuQQ9HWsyOxB7eB/cUYT5/nKzEZllLAKfNX1jJY2FR0P6p17QBEAq\nvJFvdJSZ0Sy6viUnauFWEQ6vJd7+alPz9ZP5NRZYVcifjBfMyGWqN97bcRFOdBR1mPU96gJR\nHmAQRVnWF6qmUi+o5cfWlAlsI+eXl11g+YlauKkDJWDdZUctzbeRRIGE3dq4JoliCvtSzzqx\n5GAE2QDs/ljU7y00rRAhpm1Cg64Z9YUa5LDeortKUJ59JogQje/Y2Q0XWH6iFiqJcEXVBd5J\na+AER4GokCxPiOKdQdPjxFrM2NJcOdwhiq/g/ljQ720By5QIs9zhMdK7QKadUUhBBMQaI/wt\nokIs3gPkuhCT6PDWWpqmRDh7EebkahHGy2Ha7yR9K5nqj5jmyxcgUO0ilSzAkeo4TMoNPp2j\nHkS4tEJRdTtTZGkZP8NXsIRszgKbvZEAuCfBxYKwL9QEqQm55i08NCTC2dWaiCcXmBX3WiJC\niWX4WzwKMYqUiUNHM0tpj9JYNCiC98nryp1pX4S3LtPsyhBvZYrw9LPZZB5t/bGMVak+ncY+\nkgjYMGWxc1wgK69v8YlaSEU4128inrxXKO2UJg9o+3p71h6HvENS9b/YC/IIVbLQBL+YQYAI\n113TvXE/EWaZJjicSjEEDSc1nkSAbzu//eEvbwpRMWlcxyGNB/Waq43M4hO1MBVVZTvOnNRJ\nWNFp7AozFIvBm0zC2Jzqj0lIkqBuMCPCpEc+QmHNvD+DCGfbVAWo5i7vkERQRoh8MBA+X4j6\nNd4WGo2GgjnJLddSfqIWbhVhptD1WciIMBfqgAaJrIlTVzhoLMhdQS7Nq0D+tz3mLkVJvsGw\n3BsL+p0Ncm7rncz/rBONGTWR0/8lgBle/mRDNMhJZZOlHGQl5GS4ljxzLeUnauGqCFFTBSIs\nm4bMcshpFyJN9m/YiXSlqTyREJsLZvBSZFnEY19nS3sQ4S1LdL53cF9BtzyZDcjbOC6RkyND\n8BPzanEhAFJK5RJSEfK6qszg3UUY6SPi7botLbO2hcsBQhWJaMbmRV9jKv83uj2zYaYqTHoU\nexuNa/WEdyDC2etbeOlo1og41oDwkmcX70xEIhSLp4/GiAmM1lHiEIezMmgjTiyz5JIWn6iF\njE8PcAnJXITp6yMyMonORCzDtyz9EsoEtZZBJh+MpO7FkLCi3nCqf6a63vx2IcLpVpY5SQhC\nUDr8WQwisJAX4bk43JfHknp12rp4RTsKa1GrZ/X3F6HVR1Jq0mDmxBufgTByop72YlN4E2PI\ntLMTxJZtabMszEWoJ72OknK9iTCyiYvast6KJPaQCIVjGBNFBlWhpBnCoJLBgatUICEM+8JB\ni/hQjIsnpikR6rmspmYbyopwaGv4mGguR3/qBY2OeNplOxsHi25ySvnGvV6/siIv0ZcIo2ta\nwLPYuhD5LAwX1QTqBgsfQJgf9iRgylItfpXMYkLVDoXZAhRfi7muxSdqYboDDOHKHYEaJXOE\ndWNsaFqPe0bFskGVWHf4yJtpSWBreQP5Y5+kxQov7cFEaOsuE6GIL/FZGoVYsQD9kQL/SA4b\nlo2XE9GT3OlAU4AiXLJms3Oy5EQtzIjw8lLq2dUgBrFmMBnELIENDUFnzAa8uLuJhpU7kN9k\nr24SuradmLPa1/Av3piZvbKyYvfGqn7ja4pmarZNEVQsQuVMidMCqsShI2FwOCKOzgSl0AdH\nQ2M5AlOAVyHRzfJ5aVCEMl9m/V5tJ2vNrG+COwY66+DNdL6joFEKaiaIRdR4msiUkotAY/rs\nIjRRB5u/iHGrSqWQvd54dLBpMtlGDVZfRCmDfCbu1Swfe5gNaeTXy9KH6UkpP1ELV3u+Yhej\n8tF06OugErVvZGYRHeFIMTMsKmI6ZD/bmFyyA9bmic2sHatptODCHlSEhg/S0CDWQt6E6SmV\nJB/MpSOj3uDpQmU5Ix5YJBTzFQ/cDnGFBNsWYWlDMBGRCCPPI6JAG0nBUGpcIcy0KjDI2cBS\nC0y+UB6UJorGWqrB2LDMl7o/1vULzMgyV8ej7+LBWAoBTBh4tiBsElhWKC+qYQWGeBXochBO\nlVijeLwAGeft89+gCE2wUdZUNEORERMTGZAoriemWKbcuDgyEgzMBjtAsKVqV00GEV1HzOn8\ndZWU7EyE8YrOURdETUHpA+cjCsDXkToT9eg33IwZlxgPYPygq0RalcDWuMkbrj6PJkWYepDJ\nYvlWeWqjyETsr9GqvIoDREsYMIjRyugJ1fxSxNuELSm3MCUl+xIhrvMZEYp/EhJl2uFQcoq0\ncWQjSiSC5gSwUsb/I4Xnxh9CKGawBC2KUEO7iZPDW9CpAFcItg46McbUnNV/pOfBjSF/AUgG\nh2lbD4bDOyikHxFORiykllFnj/BdykoFVoxaTTmoJ8+1MBFE4akdlQGg99Priywht1MPLYow\nZkpP4DRqkcjpEeXagGjU+EkUoVpTLhM9WxjFo7gAxCUau1CZrSyaF2FqNtGFqe0zZjUALXAI\nqrOtlOlHR0j8v4ahbGAhxtGSkXzh+iIet6C1IxEiiVFAmZbW85RMIDCitliXgtAJj6aNosQN\nVZWgUK6REva0NRoXYWo2ycyV9UnYuoQeegRcITs64xvhZgUoUO9pQOSSkT+vLRYnXCaF/Ncq\n6EmEQdVnmckUV78IpZO7A2prNQgxCmSyiblSW4vOFdYCknQfebQuwvEVzOLwLpqIWwOBomS0\njlg8sH4685ZBbpJPWF2L/NDCJi4Z3ifW5m1oUYTAUnyQNXdVhDYcnJpBtMnIofIHGSKQaJpH\nH6rDvHaN9dC2CI01jD7kFjsEmKIrUIp4Txt1ilaZC7kfEaA3qRqNRboBjpPxZ7/WQZsiBG3Y\n8iJC8DW5eQG5wLm81ZVgaPz+B/cjAmvOsKMGNGBFtRtXqaoZqbYvQgn5Lgfsh1SEIZpM5kAU\nZ7jXg2xMeSvmzxIdApQkbRU6p6hfMwY7zHuL8OXzt2aPV3u+nIxikKAzz3SoSQMfFlhW+o6d\nJSVxIJdvf2JCQYQiP14McMqOskyE6EVvZ/N6C9swWCrC2GIRnCHLIX9gsZieeNJiB4l8BGEQ\nmoYIRVJJa6S199HkBttDPL7Khq9AhKdxbcHi/HWAEQw8GyHYeZEpROMpYb8xaklJmw6MzYkF\nFcsKcZFpO5ixEBh87G/qEtUWyMXcgOv1t2GwcNzi8HC6As4vHzaRpjAe7GySMhKUGS5gFGgX\njFbX7qOFlI09I5LqcGZRIMLjv/ctWCwVoXzFeY9Ls4mLYkZK1KBCUSENdYQ/4pb4H1hO5E5t\nJKX0zVGVubZbcL3+NgwuESHJBqWqUU2RzCmRFaGZIFWVyjCgrJjBkVRxqiBcCF7ANlBMRfTN\nElkjerEoEOEZ/32+1GaxUIQ2mpiuxIGEKBFml8x0CvkhiAiFPynExhbjlyCHeGjcRBqAzox1\nwuOVu/4AABoBSURBVO6uRVn1+gzO9gvuCl0cKy0xkKINjAghDFI2LRljZfCBKHQtbXVoRIh+\nMloiBRdaA4UiPOHncBrf1x16ltMQdlwtLVZOZKgWV2hWN6mrRPmTQckCItWZjkRaJCPxGa8W\n63EXEYbqDM70SwpzlUZ2MHvoC1G+MHfs3Fi+hPNqoxhsgNjXRrYUPLO1BfwCfbcjwu/Xy4Bf\nt+95PE2E03DN7Oo/swKkG3GQoiDCuxF/3Ap0HtTZ2daRK9B4UPrs4JITUKNgJq6htHptBqeX\nDvg641bMBbPENLgQv2ltaIDSkVrRB0IzUo7HIMtCdAkBDjrBCde7McpEePw8GdGX7+OJx7fN\ne+bzxm/B0VzR8UXsGzo6bU0nf1ST7mWT1Yq+gxXAFuHT7NCyLtK699nFXICiQhswODNuWb42\nI4DZNR4IoxgmAu1cAL+JoSXhbjZzqu5O1wLUE3tsHCMUZ92B690aJSL875zWf/wMh6sNaaYh\ns9ZtgKDsYUtqJ2EusazwoOX+0AdiEIODM/ZYGZNGr09I3t3Fayw7CQWtpy1nsQmDU+1Q9J91\nZyGkV6+q0kgDDV5AYamegEH1q0ak5rtIDPUJjjGgzQVbuz0KRHhO6L+OfPiwfc9IjnE/yJ5p\nCe1bPOVs2ggTAUwigHMVrQ4y0rYZdwFHeREWNHDFTaY9TGMbBstFGFUZD6rR4zA/UUwYyxm3\nZmyo6jLbjJG1+kVWGa8m7ovFR9d4q4sCEdLb9z171vkQ6YgG8yuTDZlOO1JpD59LwVaoEqha\ngy7QHIoAp2nJSWotmeX1CkpswuC0CFUOSQAQIqs2VBBRRjoMohH9F6wCpSmWvF0BRtNqSwMv\nMFWpDEH1OXeVNVEgwuNUkS16lhkljmSMKeRCcVNRQK+zHQE30pgp5RZaluWig72Ww03Ej8Ue\nLVetjgi3YXAmlInmXA5HXIpVk3gkmGqoGjaXZBjkMxo1wTECglGBuE7GkvKNWPxreVuDAhHe\no2dY7Sg9kWMkiViEGr2M3BLOqVABd3OVDAK7jS1DgzCgeODJBaUizGvzGmqKcBtM9Bsve5hB\nDf104QfrN7GqNJS4QDgLYSU6UbDDWFJohPUhossssuW8rUETItTrZVpUSkE/iSvSeVLHFYtQ\nKRkmXLJ4oFQMItpckCGOk3Q089eTnqdVVBZb4gZFaOcKHQvPOa99Xf+sEmMcVYABGQSvCS4S\nVwshwJnyCkPFxVYCrqT2nOXRhgjlVQyjGDmcKTRu1lQZ8yoWbzwkEcyftgcGNiBrHBvBAIEU\n7XDieqgeccWWuC0RavxhZpDwHER9OmtqN4VosI5mLzuG1A3BWFhtkEMq1GpQ7oMZxVYTNokW\nRBhxxfIw2kCPxf6LP4QQNC+QxctpwR+EMLBGhKKhvp5V35fOxBxD2nUtHgvbaU2EpK9akr1P\nbFq1VmRg5YgKMFoTwCRPFXHcg/xCnEU6qdH4AhJYdZoK0JgIhaHxAPomMiJE1QQjViAZ6APy\nsQ3WIHhI9Mw6ipILguV3RyKbEqFhIbauQR2PnSYWn/FsF2gMqiSISK3DC9IsDwNPGmNrKQai\n5fWuaEOEahZNjClUapxiwkekVK3iUNe4QD3PhYPtL9ImjqtYVRANLZy+MpFP1r6h7i2YEmE0\nw2giwcmpBYwtqzZhGAxiX7UrXAZkPJ32FI/WrDC7HsRb1pukIjQhQjWNIciMBkIRKkVgMTXE\nMDSrB4x+Z73xkkHYI8sF2GcYV8F80KhAs1RKqkGfK9CECHXugAxd5MAYGDauGmDuxlLIYJDp\nsSJU48wCDOLIwJAmg0WbZ/Rtbcf90IYIhTA5F5k+tZTIqGqHNXY+ojm8GtaxHwhKjIi1UxQE\nxjElEyKDLSgOC6FY5HPd7gDkS2eSyTLhP3g7Xepw5cafBeMBg5hJ+IIzzYeBRh0RwdjmrwOs\nwjOKUCc48OKEKQxCkfot9G5Ai8Yvf+IhpUleIOORjPbAZOMA7YqB8afEgiWen7/Eul+rMNvY\n2oo3guKPoyGE9CKiShjDaqq08z9lEOIQtLi6GozkgU/VvNrugstQFd59QlsQ4fCqhhD9oPF7\nRiURuUjf2GxMIEct6mDHntmSklKo5/MizNrXTONz1/4wItRLIEOkyR3SOBQcYEgEGBlcCIqC\nkA4UB9usvl2dHSPCa15zG+wvQuOFAnzGgEZjEeSURSjs/UUUR6qW8CUeAOtQOQbzmTWPEzYT\nFF9y6cmXNWhJhBK5iAEVAkFZYncDkw0Mqg8VQwvqhdWiVjPLVOm8RpXvL8EmRIjiuhyQSJLg\nW+wA+d3QlxjEsUHpEwUino7phq7Rp+VUNcvwVSJzi2b9fDcjQqZFTR96RowzlElrQjmVTLxg\n4hDR983Sm5seS9Ct8387mhCh0sRRTNB/6CitO0T2DEPa8FBXXKqGLlLQaM/YWH1JVXWb+4qN\nr1zmTa3dG2Q/ptYyikc1excfF+uP6wTUrklDYkqDnbqY3ZCZnrjS7PzfxTPuL8IAvAVkTBUJ\nAhxLKHvyKJpOl5g28W7q24ytjEXO/QUNmKa4iS9izbWDBb9ltlsQobFkIJ3IkcG0/cUmlIuq\n7MD2KvOaMgQuMDu+9HR6fJrmG+1jIVoQIelODHg80aUNaoLVXzZWBDfGreDUw2IgkmOx+zNj\nnBj92knKjHq9EBsQoZi38X9xdRZjFaM/jDWDxDmyKuRULioFF5wbX15BS+znbSwXowERqj0D\nu8czCz7M+D81gZYDdn4aEPHC4H4Smxw5xEIR3mgjo5VzS3NNiBAmGFJB4wRjBjmwCRytjFKD\ntoDNrJxpfuay8lwgwhvjnWI0IcJgjV/i/wx5f6zYXJtiJVl7qMNkCygEXDvYSIEJnDLAq3CL\nyW1BhGq7hMRIO4bBEICLIN+swCA4sWJEK6uFS6fBRTjRgbCHwhnmFckLckICkkg+eiRK6yRQ\nQhJ1GAQtAdWbXXyEm9jeX4SQS+gkixm1BtTIiMWkxwga1iJArGzbqSQX0lVu8J5LhDqNQOKf\nZU8LqGCEfWiShJiQFBVHq4xhE+hiqzq6a+hdhHG0SBkGxxrqJVOfObYGHyhqXgPQoYDEORTm\nJiJnqssv8c5SKDmxRc9gCMk6P/khFjGzUoOEVZxYIDWAV8PAJjKb0K6xw9tevkHfIoRpHyff\nMhjHlUFJsLlBCJYxDYuAVBCchDf63Y4QzHGUcBReYrleb0ErIkwClzH6NJm7DRnVcYkcw7ga\nkBthyASi0cyCCvYQYd85oXiqEDMoEaqYRYgpKasFFakINdEtF5POs/OnVjhzsvgi7xEQNSHC\n1HSKm+PcDsOcAPrDf5pTAuGEdAVjZTODyYtweyZuMbltiDCnvkiCEupojjDVKEcs6ichkBEP\nN5QMlG+Oj+1iVRehDRHKrpmZTmI7Z0JOWK4qVSXJeknTHBjZeL0jUTGh94lJ1gt9dxEGdYEc\nh2hkCqpSR1YgwjjGBFeqe+jSdKY1odRFWNRz4r7UyMmka4BPOrEUDFeE7GMfSpca2dxoUKy5\nU02iARHCvEf3HtBcoiscW5iwhpj5GRFCy1AYKKNEcy7Cwp413pfwcziuEytWFf1YpNrUxGYF\nlHM6RpqmQPMc7i9CTBYwk9NJlTPqEKOoxrTKlIK9TfKIJObRvszw8gljU2hEhEZMYNb4NIc5\nQQkKIw8Rk/ZraSiJ7nNinM2S2IQIVYjgr+T+AerQBqYhHj86taCsqwYNH6RLJAjRBAUirbaJ\nNkQo/g0VBiHo8KLRTKLRqKnJb/ODcRGu7VcoI6EMEnwQYSKjkFwAJcILuBrytXRVwGjkUNki\n2A2tiPByhHkK4BsxRZxRXp2xzKiwXRabECHGoHxAcwOWiPKLDWQuII5nsFbWf2JHwdTrAC2J\nEFPr4R9ZY2jywS3Gkhdh49HM/iLUGBRyM+O7hFh0hldmPXaQ5kQuk8ymfz0IsTURkuUOw5E4\nHuE6FUZ6JeRsm8gGRKg6RHcHWRwmdBqcLosxtHTKh4RHUYOtG9ABLYlQItCQijDklVJrkjcL\nOe+g391FyCkYSw88Xgi8tSKq43A179Nm5muWa8oXaT2VGNCMCCFd4PAlkLxjBSPCcOs6t/Fu\nZdzFELciQhAdYTaB95/A/UW7NMPR+fkqUGi0KxcPtU00IsLIZkbZoXGF0SRDyTWDwPhoXROz\nzcPrZmhRhOMBOQ2hJ5l3O/hb5itHoItwSc/AI+6zcVwzFkrsJARABRrK0LStSu6zBnYXoc6i\n2XIxaYROBQvUCNW2WG37zUW4oGdRmeT27BDJlFRtcj2NW8F3ZinMbqpF+3SV8TQilFknoVDm\nOsklTIphuJo0tjeOsXENtiVC8IDMUYjpsIfEbYo48zvYFCeWti0X4Y39ytyj+UyliAnEUC1g\nDGO3bqqMsZ6eN0RTIiR5YhtlmasGIgTXGRWQRARZt0tnLLptPLrxVLYgQj04xZ4hR9QaMNKB\nTGSy/eWjbF6CG4jwv8+3yxS/ffy3oGfNKvB9phZaz4DuLPfIU+AScTCb3R2oiBlDXG9x1B59\nKYPZmH/yXDIVbDxRdYag9sVTC5VFeHwhxWt5z1YxdH1YIMLAXI4nwP2BnuMbHiEOjbbBhNZq\n9lt39OUMLhNhLqixISuHKnEbPfiy21BZhB90+Pdz+fT7faCPBT1zGBOlCVO1QEl2Nya9qyix\nrmoc2qiwE7eihZpxat0lWs5gmQindslGIiAUgRQDvvaR1d2GyiI80I98/qHDop51V22u73j5\ngpJs2iHnUYR2tzWVwho9rVglcw5jMeou0HIG8y7enspPjt57YvUpLdHe6n3S6p1RWYSZfQ97\nGpCpmTme7SHd19ZwNhIyxqN0ZYmssrprVknDIpxncIZAqJHfjoEmxoAn2p7WOCgTk66+oA7Q\njCdkj3a9j6iQ2QhHHUqzydbNVFt301PDIrzRE9opzUWnGrYQK1KjUgKeJlp4QFQW4Smj+P69\nfFqYE4pOVoT/U2kHv0iJjKebXzIlna+p1HZOWMhgib1MisKDFRiD8tk4VnERrsErhCsvx/Ke\nSf7Vm27d7KaYZfkaQJX3FOGawPfKCCqhmMGyoEUd3fiB2AmaKAWjl4wrfGgNVhdh+O/jcpfp\n8PZZfp8Q9lWo/oTP3awzIzFes7x508aCUdW6ztrzVcrg9X4h+9M6okOwkCLFtGnfHb1Lz8gA\nmsOyVq7zM62SjAu8cockO4Kw7yrZq+sCESYZxhjswH3d4QufzHHl9wnv0LM4ILIxScH0lwhg\nJl5MRJjZGSjAzqukWRGOeXi89yV7M8PpoAzS/hZtF+wvQjWIeLcvJAlbjpqSUHCBCLWXnpZB\n2yKMCsZ3gkkOjl8e3++laEiEcSiIcpi86UvQRFE/mXPJNoCLsEq/+ck0KmMC+jN9NdGKCCkx\ng+Z8lqTo/t98R9kyE3f++1oOzYqwQFoYlHY16VWxvwjTiMR+ojClDfPQ00xHc3lGHP10aJIb\nFuHVDC+XGT4fWhChzc3j85MiZPUV7I+W89vhcmhXhAUzL3FQX3NeFw2IMNhd6qiACVQSERr9\n1uGxu+XQsgivN9Kf0auPJkR4+ZqwYZ47y1QZ3aA8mRbXfxZ0LcIOjV59NCPCHBv4uFlIJWaU\n2WEyVwnNiNDltBINifBK8ZTh1FWCZ3weNCLC5w1FbkZTIlwsoOTnzp5xJbQiwtxBRwkaEuFN\nArI3nJ7KHbYhwvxdJEcJWhLhTf3C3Qp+eP9J0IgIO3zgrxW0I8JVpjS7dWN/m+XDowURqt17\nmmmviK5FmDzizTcVzc3FR0cTIoT/HUvRtwiTCoMtfq7AqAERsjF8qly8HtoR4XLpZGSLK+FZ\n1kMjInQJrkZLIly6O5r1nfRsGmxHhE8053XRkAgXm9Ip6p/sZmEDInyyBKA2mhLhuiZq3PXv\nGU2I8LnsXmX0LUKnPrQhwieze5XRtQid+jPaEKFjPToXocNF2D96FOEW7q9jl9qQCDuexT3R\nnwi3SAS7Ti6bEWHXs7gnWhTh5F/dgqfRKotwgzbvhnZEmD/suIb2RDhhT/nwFveF+77X3IoI\n+57FPdGgCPNn+fCeImwz5XER9o7mRDhBpRzeT4Stpjwuwt7RnQj3ywlbTXlaEWGzE9Q8GhHh\n1T9SDiLcaXe0WUPfjggbDRWaRxMizPwO7jkzu899Qhfh9X7bTJqbRxsixANXdkd3g4uwkX4f\nDy2IMIlEZ+8T7oZWUx4XYe9oUYSNYndfPAEXYe9wES7A3r44Dxdh72hBhHOBXpvrvim4CHtH\nGyKcCvRajQCbgouwdzQhwum9mPsMpW+4CHtHIyKcL+J0z8FF2DtchN3DRdg7XITdw0XYO5oW\noeeEJXAR9o7GRei7o9fhIuwdbYvQ7xMWwEXYO1oXoeMqXIS9w0XYPVyEvcNF2D1chL3DRdg9\nXIS9w0XYPVyEvcNF2D1chL3DRdg9XIS9w0XYPVyEvcNF2D1chL3DRdg9XIS9o3kR+nNr1+Ai\n7B2Ni9Cf4L4OF2HvaF2E5UWfFi7C3tG2CP2negvgIuwdLsLu4SLsHS7C7uEi7B1ti9BzwgK4\nCHtH6yL03dGrcBH2jsZF6PcJr8NF2DuaF6HjGlyEvcNF2D1chL3DRdg9XIS9w0XYPVyEvcNF\n2D1chL3DRdg9XIS9w0XYPVyEvcNF2D1chL3DRdg9XIS9w0XYPVyEvcNF2D1chL3DRdg9XIS9\nw0XYPVyEvcNF2D1chL1jGxF+vRC9fa/s2bEMm0xkAYNOYC1UFuHww3+vdMHHup4dy1B3IssZ\ndAJrYQsRftDHMYTfD/pa1bNjGTYQYRGDTmAtbCHCAx3Pn4/0sqpnxzJsIMIiBp3AWthChPwL\nKeZ/MYVzWAkbiLCIQSewFrYQ4TtTeFjVs2MZNhBhEYNOYC1UF+Hb59c3/Tt9PH54Xn8X1BZh\nKYNOYC1UF+GAy8fDcVXPjmWoLcJSBp3AWqh9n/Dn5+vr7e2S3H+kDBJiXQeOGJUncpZBJ3AL\n+BMz3cOfmOkdLsLu4SLsHVuJ8Hqw8hActhCUbTSCq1e2/5U/ClyEN6CNv5ThIuwdLsIbQPC6\n9yjqN+sivBdchOtB0ftOcBH2DhfhergIHVXgIlwPF6GjCvwWxQ146Jyw2X4fDy7CG/DQu6PN\n9vt42FGEjkrYmikncGtMTvFd+fTKd6y8BTqdisYruwgftvIW6HQqGq/sInzYylug06lovLKL\n8GErb4FOp6Lxyi7Ch628BTqdisYruwgftvIW6HQqGq/sInzYylug06lovLKL8GErb4FOp6Lx\nyi7Ch628BTqdisYruwgftvIW6HQqGq/sInzYylug06lovLKL8GErb4FOp6Lxys3R7HA8G1yE\nDsfOcBE6HDvDRehw7AwXocOxM1yEDsfOcBE6HDvDRehw7AwXocOxM1yEDsfOcBE6HDvDRehw\n7AwXocOxM1yEDsfOcBE6HDvDRehw7AwXocOxM7YX4ceBDh/HuQNLKn+93FD5hP9KLzip/PNO\n9P67rvJxyTWfLtIOclnl+nAGN2ZwcxG+Xv4o1MvMgSWVPy4HDkWTkevoeCi84KTy9w09/x6G\nyoUL4Mf+Ea0lE7YFnMGtGdxahP/R4Sf8HOi/yQNLKv/Q+/FsZt7XVD7jrfDP/KWVD6cDxzf6\nWFP5/VLto2jY4VwTB7lkwraAMxi2ZnBrEX7Q9+n1H31OHlhS+W0YbhEPuY7+zfylxvnK/y4s\nHOmwpjItGPZpgb6agksmbAs4g2FrBrcW4RudPfgPvU0eWFJ5RNFcZCr/RrOzoPI7/RRVzFYe\nI6gi/sNprZhBLpmwLeAMhq0Z3FqEiQlZYlMmyh7pdV3lV/otpDCp/ELh83AJpVZU/hyDmSJn\n9hNd8CIjvAGcwbA1gz2K8Ovi3pdX/qR/pWs5M+y3S2a+qnL4Ouf1h6+irsMTiNAZtKVLm12J\nDSj8PRRFZknlSzywnsJzWv9eZApzi+eM4qzu0UXoDEalS5tdifoUHg8loUwuHjnvTq+n8JxR\n/BZtzCeVv87BzIn/UkP64CJ0BuPSha2uxSEeQHJgSeUzXgvvmMWV3y8hUCGFSc9LVl5S+YXO\nmcix+F6f6WXJhG0BZzBszeDW1A47Q7/x3trvgr01U/b35bXwhmlcmQTrhj02sqbyUmdmCi6Z\nsC3gDIatGdxahJ8X6/Wtt0iTA0sqnz6XRTKZyosonBj2b1H3SeXBFJbdohqGOtfaneEMhq0Z\n3FqEdZ+3KJvCuY4KbVmm55fjOSn4t6byB50fHPwo1pEZ5GM9MeMM5koXtroaLxfT9aojgwOL\nK78vMIVpz/bTwsqfNwx7fHiwePnxIJdP2BZwBsPGDG4uwuHxcxgRHFhceUk8kvZsPy2t/P26\netjjY/RllUNM4ZIJ2wLO4BlbMrjXnpvD4RjhInQ4doaL0OHYGS5Ch2NnuAgdjp3hInQ4doaL\n0OHYGS5Ch2NnuAgdjp3hInQ4doaL0OHYGS5Ch2NnuAgdjp3hInQ4doaL0OHYGS5Ch2NnuAgd\njp3hInQ4doaL0OHYGS5Ch2NnuAgdjp3hInQ4doaL0OHYGS5Ch2NnuAgdjp3hInQ4doaL0OHY\nGS5Ch2NnuAgdjp3hInQ4doaL0OHYGS5Ch2NnuAgdjp3hInQ4doaL0OHYGS5Ch2NnuAgdjp3x\n7CJ8pf9Or//R+94DcazEAzD47CL8pcPp9XA47j0Qx0o8AIPPLsLwRZ/hk/7tPQzHavTP4NOL\n8BTOfNHb3oNw3IDuGXQR/hLR796DcNyA7hl0EYYP+th7CI6b0DuDLsLu7ejTo3sGXYRvp4zi\nde9BOG5A9ww+vQj/nUKZT/raexiO1eifwWcX4fFwucvUdzjz1HgABp9dhO/j8xZdhzNPjQdg\n8NlF6HDsDhehw7EzXIQOx85wETocO+N/dANFk+EnEFgAAAAASUVORK5CYII=",
      "text/plain": [
       "Plot with title \"normal errors\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "options(repr.plot.width=7.5,repr.plot.height=3.5)\n",
    "par(mar=c(4,4,1,1),mfcol=c(1,2))\n",
    "ylim <- c(-5,10)\n",
    "plot(X[,2],Yt,xlab=\"x\",ylab=\"y\",main=\"t errors\",ylim=ylim)\n",
    "points(X[,2],X%*%beta,type='l',col='red',lwd=2)\n",
    "plot(X[,2],Yn,xlab=\"x\",ylab=\"y\",main=\"normal errors\",ylim=ylim)\n",
    "points(X[,2],X%*%beta,type='l',col='red',lwd=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The t distribution produces a lot of outliers (some outliers are not shown because we restricted the y-axis to ylim in each plot). Given the data, our goal is to infer the parameters $\\beta$. We will consider two methods. The least squares method we discussed in class is\n",
    "$$ \\beta_{LS} = \\text{argmin}_{\\beta} (y_i - \\beta_1 - \\beta_2 x_i - \\beta_3 x_i^2)^2 =  \\text{argmin}_{\\beta} (Y - X\\beta)^T(Y - X\\beta) = (X^TX)^{-1}X^TY$$\n",
    "A second possibility is to find the [least absolute deviations](https://en.wikipedia.org/wiki/Least_absolute_deviations) estimator. We call this $\\beta_{LAD}$ which has the form\n",
    "$$\\beta_{LAD} = \\text{argmin}_{\\beta} |y_i - \\beta_1 - \\beta_2 x_i - \\beta_3 x_i^2|.$$\n",
    "This may also be known as $L_1$ regression or quantile regression with the $0.5$ quantile.\n",
    "\n",
    "In this homework you will fit both $\\beta_{LS}$ and $\\beta_{LAD}$ and compare the methods based on parameter estimation accuracy and speed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 1\n",
    "\n",
    "Fixing $X$, generate the $Yt$ and $Yn$ data $N=100$ times (note that $n$ is the sample size and $N$ is the number of times we run the simulation). Compute $\\beta_{LS}$ and $\\beta_{LAD}$ each time. Since $\\beta$ is a vector of length $3$, this will result in 4 matrices (LS for $Yt$ and $Yn$ and LAD for $Yt$ and $Yn$) with $N$ rows and 3 columns (you could store this output in other formats as well, but matrices are convenient).\n",
    "\n",
    "Make scatterplots of $\\beta_2$ versus $\\beta_3$ estimates for each type of data, estimator. This will produce four scatterplots. Mark the true $\\beta_2$ and $\\beta_3$ on each plot and comment of the quality of the estimators in the different settings. You scatterplot should look similar to the ones below.\n",
    "\n",
    "**Tips for Fitting** For fitting $\\beta_{LS}$, you have many options in R (lm, lm.fit) and python (scipy.linalg.lstsq, sklearn.linear_model.LinearRegression). For fitting $\\beta_{LAD}$ in R you can use l1fit from the package L1pack or  rq.fit from quantreg. For fitting $\\beta_{LAD}$ in python the options are somewhat more limited:\n",
    "* [statsmodels](http://www.statsmodels.org/dev/examples/notebooks/generated/quantile_regression.html) requires use of pandas, calls the quantreg R package\n",
    "* generic optimizer such as [scipy.optimize](https://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html)\n",
    "\n",
    "The later is a more general method that will likely be slower. But scipy.optimize is general and will be useful for many other problems. You may also use any other optimizer you find. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Installing package into 'C:/Users/Kanchan/Documents/R/win-library/3.4'\n",
      "(as 'lib' is unspecified)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "package 'matlib' successfully unpacked and MD5 sums checked\n",
      "\n",
      "The downloaded binary packages are in\n",
      "\tC:\\Users\\Kanchan\\AppData\\Local\\Temp\\Rtmp04U1sr\\downloaded_packages\n"
     ]
    }
   ],
   "source": [
    "install.packages('matlib', repos = 'https://cran.r-project.org')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Installing package into 'C:/Users/Kanchan/Documents/R/win-library/3.4'\n",
      "(as 'lib' is unspecified)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "package 'L1pack' successfully unpacked and MD5 sums checked\n",
      "\n",
      "The downloaded binary packages are in\n",
      "\tC:\\Users\\Kanchan\\AppData\\Local\\Temp\\Rtmp04U1sr\\downloaded_packages\n"
     ]
    }
   ],
   "source": [
    "install.packages('L1pack', repos = 'https://cran.r-project.org')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "\"package 'matlib' was built under R version 3.4.3\""
     ]
    }
   ],
   "source": [
    "library(matlib)\n",
    "library(L1pack)\n",
    "N = 100\n",
    "Xtrans = t(X)\n",
    "M = Xtrans%*%X\n",
    "X1 = cbind(x,x^2) ##required for l1fit\n",
    "Y_t = matrix(0, nrow = n, ncol = N)\n",
    "Y_n = matrix(0, nrow = n, ncol = N)\n",
    "ls_yt = matrix(0, nrow = N, ncol = 3)\n",
    "ls_yn = matrix(0, nrow = N, ncol = 3)\n",
    "lad_yt = matrix(0, nrow = N, ncol = 3)\n",
    "lad_yn = matrix(0, nrow = N, ncol = 3)\n",
    "##Computing least square estimates by directly using the formula for beta = (inv(X'X)*X'*Y)\n",
    "##Computing least absolute deviations by using l1fit from package L1pack\n",
    "for (i in 1:N)\n",
    "    {\n",
    "    Y_t[, i] = X%*%beta + rt(n,df=2)\n",
    "    ls_yt[i, ] = (inv(M))%*%Xtrans%*%Y_t[, i]\n",
    "    Y_n[, i] = X%*%beta + rnorm(n)\n",
    "    ls_yn[i, ] = (inv(M))%*%Xtrans%*%Y_n[, i]\n",
    "    lad1 = l1fit(X1, Y_t[, i])\n",
    "    lad_yt[i, ] = lad1$coefficients\n",
    "    lad2 = l1fit(X1, Y_n[, i])\n",
    "    lad_yn[i, ] = lad2$coefficients\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4QAAAGkCAMAAABeoxsWAAAAM1BMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD/AAD///89ODILAAAACXBIWXMAABJ0\nAAASdAHeZh94AAAgAElEQVR4nO2dh7qiOhRGg/04ys37P+2VElIILaTt5F/fjEdFSGSzSAUZ\nBwAkhaXOAAC1AwkBSAwkBCAxkBCAxBQgIWPGl3jfG8Yuj8/6avdwOQK7mAWO88fvvYeyvOf6\n+K5viHooC5TwIYL3XFnpX1PANyeORcI+btqLgffKZuiHknr++SyWLxm6f7tXAgmYx+CvD9uf\nXF5HKKnnn8+CcPlVaFrOv1fGbrtXAgmYx+AXs1895qovb9+/kDZHNkMN6vnnsyCIl+1KcMbT\nq/JO+2hYM7Y9fku+l65pMj3p2pm/F/c31z/QPn+HDbv9mdsHe5gF6NvJ9jPuay6/LFdISwhl\neRL+AnZbqbyIVfTIfRul1vP7e+nPx9OT7hTdc+PqB8Ra07kbHGEm4aNrxz+nrhm5/P2zZnkj\n5ENZnoR9m7C5/631qM0iJ0LQyKV/ypPb1Da5qR+49wvbX1hfAb5X8cwk/EWh7Wowjbn8995l\neSPkQ1mehNNuvqx1qekr/cS9/sI/RKJb9nulPPmdiNnrt/x3ku6rReJ9NtScVg4RsIIZuPdQ\nDt1E3ZNpHaU7N0MxlAVKyLuGfM/DvsJ8pVt3Cu7fvQ3LpqOgf3IX58fHUC0S73cn3fua6mAN\ny9mz25lvpaq49NHlzVAMZYkS/toFf/e+VvJSP6TVWvSVpipKX4lhYxynJ9Mb32E18fo5rAMP\n3TACN9VDG7nbbR8tL5RlStjxvalVi72RY+oy8UR+2PiAmBfQbEzpADaMwCkDvMPZUy5vyw5l\ncRKK0yg/cvpszDjqT9TTZ6Ot3P4NnW3oHnXAkFA0Ijou+nKtd7S8UBYn4X2K13ehiLSsdNPG\noWaRu80bEsrK/biTc+4rRt9t/5jKP235dW3iGv1QFnD0GPuxq9l3VYp3szy4NKzUTq/+fqfF\nf/2fK7dEztKl1r9/mRoaKxM6wBJ64B5yvtpr6FMTy//dNmfM0A5lGRKqzYCb8nK5ft/12zz0\nl/Mz8PTkOm3xrr7f9al/uTK8DI6gB47JdsQ42UldvtJhQj+U5Ukod3OzErmu3qHU/t/jKo9x\ng9OWx+VXNXDz1jyahC5ogftT6y03ZYB9y8ECQlmghMP1hOz2bNfWuumV1fbRTXczBoll5IZt\nKhMOp7e7uGG+jBNa4LRm31tMNRu8eKxGkn4oC5AQANpAQgASE1xCBjwROlIIYGgWd3HwGIZO\noBaSSZgo3fKAhOSBhNSBhOSBhNTxLuG/5zAofntsXMWOGHoCElLHs4StOtl2fXgTMfQEJKSO\nZwkfrPkb7q37fTfrE30QQ09AQup4lrBh8vbWn/Upr9YLS8BxcpAQATyDZwmNa7b2pTxOynVK\nEGQgIQJ4jhxKQnYiQZCDhEnzQR//bcL3cJnQ/jYhM98Ah0guIQJ4Et9DFFeld/SyOscdMfQE\nJKSO/3HCRz9O2Nyee8cJEcNzQELq5DBjBk2KUySXEAE8SVwJ7VPH0bl2igwkPBBADGbMyaEk\nRGRO4XnXtY+uU/vZ/RrK+i8RuQQQ51sbeUgITuB3R3a/QsTbZsfEQ5d0UXG1AQnJ43dH3tmt\n/T10N4f83lcHmRzSRReOFe8zZnZdMeyeAJjhd0f29xccbzLYrk63gIS+8CzhCxJGx7eEvJv4\npLzQl+6M7sLGjb+gx3d19NPsvT8jAuEJ39XRT3er2372YbvaKESb0Bfe24SfvXcqRiQ84XdH\nfljz+PBb8+l/znH151Qdto7eURv+O2Zeyhxup5TBMTzvyHcjK5xP/+liNGoOekfJ431H/t37\n+yPcnuu/yOeeLkTUgYTkST9j5uB6qJIaQELykJPw1NolAgnJQ0xCDFPMgITkgYTUgYTkgYTU\ngYTkISYh2oQzICF5yEmI3lEDSEgeahJinNAEEpKHnoRABxKSBxJSBxKSBxJSBxKSBxJSBxKS\nBxJSBxKSBxJSBxKSBxJSBxKSBxJSBxKSBxJSBxKSBxJSBxKSBxJSBxKSBxJSBxKSBxJSBxKS\nBxJSBxKSBxJSBxKSBxJSBxKSBxJSBxKSBxJSBxKSBxJSBxKSBxJSBxKSBxJSBxKSBxJSBxKS\nBxJSBxKSpwoJi75VKSQkTwUSFn7TbkhInhokjJ1gXCAhecqXkBl/SwMSkqcQCVdafZAwesrg\nGEVIuNrqg4TRUwbHKEPC1U2iTRg7ZXCMEiTcKOvQOxo7ZXCMCiTEOGHslMExqpCwaCAheUqQ\ncLXVV3Qp2AEJyVOGhIutvsLbgx2QkDxFSLhc3hXeM9rhW8LvnTVPzl8X1jx2p1x8fSMohUi4\nkUrJR4hnCdvmJxR7PbtHdt2Xsnt9A/J2QELqeJbwwX7l36Nh95a3/fMdKbvWNypoLOyCqoT7\nTqGQ8DBNvxpjbf+n2ZOy816uoLGwC5oS7j6FVhBmzxIyJh8tu5ipGAkdTbCGU+QuiEq4dxsV\nVHgClYTdYxu2JISEIyQlPBK94pv+gdqEj3Z8viNlx/oGJBwpXsLiIdw7WkFjYReQkDqExwkr\naCzsgqSEOIUqkJ4xU3xjYRdEJTx2Ci061KQlBB00JTzkVeGVHkhIHqoSHk6q2CMGEpKnfAlL\n78WBhOQpWcKhygoJo6cMjlGuhKIpCAmjpwyOUbCE4hFtwtgpg2MUK6EsANE7GjtlcIwKJMQ4\nYeyUwTGqkLBoICF5ipWw9KbgBCQkT8ESlt0UnCAtYdENhd2UK2EtESYsYS3nyS2C7IHXhbHb\nO366VUJZQj+bIY/fPTCc1a7DHUj2XZXtOX23y9ooHwd0Jaym72yLABI++lsjfB/sFS1dJfU1\noay2Ua8TUZZQn1hYLQEkbIbb5bXsEi1ddZvLB2UfdKa8ZjvWyh7CEhoTC6slgITiOJ+VLrbb\n5flM3Pg7z9wsdyVMLSUuIeVKiC8CSHgXh/na/fLiS8j0Mo9pb0JCnynvXZ16S8AXviW8PV9v\n9vd7unEP9QQSagU0M96FhB5T3r86Q22U+5dwqmsy1rTR0lW3uVwbVf9K99Am9J7ysdXp7npf\neN4Dn8/rdbv13TOPNQdT9I5qFdDJSEa9TkRWQvKnP2+k2gUJxgkZUz8xvgjWTRQNwhJaTn+0\nY+FIWRJupNm3QJRaKflSsIOuhHPlSoiHAzVJaAR9KgFpB52yhPYt0o6HA3VJqLM4oEmKgiSs\ntasGEkJC7ymf3SDpeDhQs4SzAXuSQELy1C2h1l9KlIIkRJuwjnR1tP5SohQlIXpHa0jXhHo5\nWJaEJcTDgdolpE9ZElYJJJwgehKGhOSBhCNkmyOQkDyQkPoF9pCQPLVLKO55QXeIqlAJiTYO\nnICE6lghxcgXKSHZxoETdUmoSCZu9CXs4+qFFZQoU8LQCWRFTRIqp9fpNk/T9FHW52n4R0rF\nEiVUo1IBVUkoH6dZozLcQznIZQ3VWDnXA6JYCUuYzrSPiiRU+l7kU2UOt2wYah+X7+V5QBQq\noTgpZrnPPVO9hKLskyWdraM04zZKiRKOgWFqrAqmdglFpcdoLqqf53Yvc6FMCZnqYelUJKGt\nTdg/ikqPcWUTJHRN2ce2md59XTRVSTjrHTXvMqP+QInRb6p9LidKlXAoCvNsh3umJgm1Lk7t\njqO6ZLahDGVxZpQrofEDPuVSl4Q60qxZn4x6d1Kl/yZm7vZSpIRMzJyowsKKJVTNWyjpZKMR\n44RRElDbBwXcHX0ndUmo33dU+bsQbfFutgbyoiSctQ/kBS45R+A0NUlomGb0xtiizKb/+R4C\nviVsH91P2j0vjF3/HFN2xey+1iYVFlwiViWhkbLW9bK0ivQwUzxL+G1+O6NthubY1S1lR+xV\nE6UPu1QNK5JwNs6g9oIurSNnlIbNnDueJbyzW/t7uH9/Pt5Xf2IyqISca7N5uRy0KJCaJVTP\nrcuHcu4/WOFZQsba8eFXM139seXAEg7PZEuBiX/lUbeEC29aVss3+N4l/D00WseIvljBKYG1\ntJVHmf4woyn7ZsEJKpLQqtOahEyrC+WK9+roh/Nn99CVhKuNwoC9o/IdZbrEasOBMlVJaNPJ\nYqb2o9889w4BzxJ+WPP48Fvzs/B9YW+nlN2Z/WCheHtqD+YcCldqktCq08xMMVuKSsx9D1G8\nG1nffLqlfICtE5w4R2bfND9DXRJaMY6DsQNg13zRHMpI3xJy/ne/dAbenl/HlHezXdWfaqMl\njxVCQgOl7bFVFOZxWPiX8HTKB7ewURayaepu8n0dCEhosCShrSbL1Q+kgrCE213T8gPlKggJ\nZ0zyyV5xbi/0dh9CYdkp4athl1eklI9uYLW2wTjXrmIpUMZ9X8h/BPPdkWObUPfu6OhGRLYk\n/NxY8+LPHdPQvKV8dAPLX0F0To/7v9CW4cbXCRbBfHej7B2VNSDr0UJDwk//NR7s3vLvjXk9\nk4ZvE4qpa5OD2r0PimH964SLYM67cdTPMrXUVhSm/iYbEvbzPx/D/LOWXaKkvH8L6+Wa7Bod\n02M7ussIsv5twkWQwF5UegMWJMyiarQh4fgFbsqL4Ckf2cZalqR7cqCCq031QtjqHu4fA0Qw\n/71oKQrn/aPJFdwp4d9Qi1mfkO0tZV/bl71iYg5TmfXRPRKGiGD+e1Ebrc+j0LOyWR29t+Mb\n7calSd5S9rb9KQJT01CpkpbDVnU0VATz343KeXj4k2mWNyTsrs8Vb/gtCGNIqI1MMPGHsXzD\n4cL6VwkXQQK7cGqRcMol4a9NLwLXeC0HI8RQaZVzcRo0+8wsa+UZpxW28hsqggT2k3IRTSYd\noVY2JYyfsrcENNdkjXQtHBmfLhdJlVsKe4nNOuNyzHXBEhqlmmwhjotsGcj4dLkIJFxGuaJQ\nvJMwN0vslPDflTWPduGjnlMOxBQPtlzg5RypRfZl1n8EaewkMWsm59BuSfi5dpOehmkXDc0Y\nyplL04CRcsmnVlzmHKlFNjIbLIJEdpI2Spxnpjck/DdMero2H95eaQ1RiGTmhZ520wvGuNZw\njJk3P6xnNlwEieyk6VTLptnDubEh4XWY9NTfqGLj7mneUvaLpZUnJ3YLBdnKp7NnPbfhIkhj\nLynVnWzn729IKEbZlBfBU/aKrWyTl1NIEadFPMswrbGe23ARpLGXlBMuV8cqchqLqlTCsZWg\nDSQpC0kBCdcYtBunLMpZpFmdbauUUMxjkzPb6KkngYSrMDlxUZx/eWbtjtIltOxt0SHDpiJw\nqpxmEpRjQEIr6nWkQsLBQLX5kS5/CpsSakRJ2S/zeocyb43J0yOffUUqTm5JGCqCWe8dNexT\nkzDXofviJZzJNHWUKTUUJt+XH+Pcb8kRCkhoQasAjZUBqhKmSDlYglrDQJm9JvpINe2yajSs\nkiqPOe8bI4pKV8w0PqV+LDHVSKjsfKXMU1/LrjNuRjFrIOEctQ+mf610wik9NJl8hXokVBNl\narVTXHE/r8HwvI+0EUg4Ry/wxHtc9ATw2UL7O5GoRcLJvlniU6NdawVCwmzT3YW1ujmG+L+e\n+ULOE5WN1Us41EpEr4U8GWbVaFgFEs6xlYTi9X8jpoTKY2Qqk3C+q0X1ZOo+VIfUcmk0rAIJ\n55htwv69IaL/TZhraH9jUouEWkNA7awX7fPpp5u0xjwFQuVy6+vnvHdsSg3P//vPbiEkjIBs\nlstu6tE52TeTJGdngYQW5pXLuYOQMNE44ZT0NC4he0chIZfV8u3B/ax3k22iVMd/FguVcWO0\nCWMhKp36LxeiJOz51xQhoaVXpn+cSyh0Re9oVNj0f6qGypMhlZagxHN+2xu79j+ybNsRoabA\nxUAZntAklAsxThgP6dsgotZDyskdXt5z+8fYHyfeJrTQh3smYQ7VnxolnLyTPaNG72m6rDng\nP7ffK7u1xUnYR3ehIISE0ZEXT8gbHuijE8my5kCIzD5Z8y5PQm4UhWqtBxJGR05SY6p/OUTk\nMEEy+7ls1srp7CRtqrB0ULvmMFXeVhMvW8I+A9NMNXP2Nq12YaCc3kuR0Oj2ZLIctC1OQb0S\nKnNnpHPZxOUAqfJJa/8ouZ3qouPy5GfciiUcH7V7wmZTQzkAJFxl3sRQYy1eJhWxWglVB9UY\nZNJWPwAkXMUSTfX6US77yLUvFNNKSDiN2xuLEp8d9wMJV7FJOLZEpvdnVZ+4DRJIyNS/YpG8\nvjBN7o4ACdeZNy70m6/Zqj5xGyTVSqgUhUqBOCopz4Opc7kDSLjOrFSTzsmRKm1B7KGqiiWc\nzoPyUorpHb22kjeQcAujXaEqptV3IGF8lKuZlEF7Ng3oEikKIeFB1JKQKxXSWf+p+RVDdRLU\nLGGP7B2VrUTb3MJ8gYRHmeo4bPwphIUaq/4Nw3XW1C7hWPqpZ0PtLiR55HIVSHgAWf3U5g6b\nA/Y24cI1T4JJuHnKyCOGojUge0S1Kb7oHc0v3ROo9zjhXGuCzD+qjmPwoO1ESDj8kdFZvAuJ\nulZOckLC3RiFmVYo2j5u6zrNX0Km45RyXMbJo+NTbrv2mhvWZTa1FBLuZeaRctPn5c+Tk3D/\nHUoyiSFTKyd9ng0JJ0fNikke+eeQcD9zCadpizs+TqdNuHqHEr7b0GiMbcH+KdO6RuUwhTFi\nGPCU6AQk3Iu9JBS9o1sfp9Q7uvMOJXnEUN3PY2+ZISHjZqsAEqZN9wRmm9DsfLF9evb5UNk6\nsmCTfXcoySOG034ez4dsdntYxnnmP9gECXdjtiuY+sfyceUxJEF6R/fcoSSPGEqh5Nw1o1+G\nmRKiTZg23VPoPWzTX3sZF6sHLoiEe+5QkkkMpVCiZqLdhYTLfmyun0AzadNySOjOVBIuBjRO\n10UYCXfcoSSTGJoXTGgWjpVRc4w3Vmx2Agld0duE6b5PKAlPpByZSSjlIgpZFx07ZrTZTpkB\nCZ2ZXWKfKhtHFwRPOR1sXj9Rx1IyLAU7IKET8wFgSJgPmoHq++PCpWGlVEBCB6bpaigJs2Sp\nhSAajiyzvENCByyDTgvfJ0LNBxKaLJ4XmfE/FyDhcYzRJrOxL8WL0g8ACU2WJGRMKQVzyjwk\nPI4ZZG0epSpelI5TSGhil1BccGiOGWYAJDzOXEIuCzxFvDjNRUg4w3ry03tGs8o7JHTAiOOs\n2DOucIKEkbHN6lWqLZyjdzRpul7Qg6zWSvXZGVz/QKDMHF4QPOX0zFVTgoFxwsTpekINsnBR\nc08pCtEmTMBs18c5I7oBCY/B5n0uqoRi+dgNh97RRGi1E/WtWZ5zKBYh4RG00m4aqp/agcz8\niBbiQOGGhBbM2glfOCPm0UCEhEeYzqbapDVNPu3+l5qCPEy4IaEFrStU75PRwpJHVykkPICs\n2ZiXastaqnofaM7l5TPBwg0JbTDRUcaUF+NZksvSUX46JZDwAGoP2/hE9VIsEtEe3gg9zxsS\n2pimaYuGgwjEqCMkTJjuKVQJNdPUgk/ciMzeIQAJY6E1C+SJUwscJCSIFGr6hUJliX59k9o3\nIItCSBgLxmWktJfDQvXcmPqLQMIjzEYHlY4X+QE12rbpbJ7zdHhB8JSzQGmRKwIyPSzoHaWH\nMlXbDN+86WfUQ0OFGxIuMNY+p3ro+M84G2KckBaGR2b41LtSy0qpLBsxThgXo9oibzbDw5wN\n3YGE+9mqUTL5K2missNm6gbK1ZEFwVPOhfkdoMTZMK+sQ8LdbHalybviz4rICNnavyB4yhmy\n2IzIAUi4my0JxSUUkW8jBAkPkls5CAkPsCmhUgeNCCQkDyTcz0ab0Oj9jgUkJA8k3M9Wc0JI\nCglzJb+aaA8kPMJ6ELUxiWhAwr3k2CfTAwk9otzMK2KihxcETzlTZs2JXH5uGBJ6Rb3naKSv\nCAl3YnasjafMDDSEhEGQ186ET+rwguAp58lMQjG/O1F+JJAwCEy5mjewi5BwJ4aE8srC5N8D\nEoZATJpR7zsTMK2DC4KnnCl6m1D+vmTy7wEJQ6Dc42L8a/uQHzch4V6MC12mGwQl6dRWgYQh\n0O9twW1f11sRCQn3o97qSbUv8eAFJAyCuJZiRcKF948ndXiBL0jGUN6aUukd9RYMxzxVlm4k\npovXxOvZB5YWHE/q8AJfEIyhfo2h/XZACXJVWbrR0HpGlwpCSBgZdVCCKe9w5Z341JZuRFZb\nG5AwBVMR2P+dQgMJS2eh3w1twgToEpr9ZmgTVoTfDjlIuB914EjNPXpHK2MKOMYJ42PUQ5UF\nGCesCa1XwEPoIeEe5NC8HKXIBkgYGaUTwE8lCBJuo45M5OVfDySMjJTQ+OGusxvcv8AXdGKY\nuONlC0gYmUlC2TvjZ4P7F/iCTAxTD0FsESZfr4ZdXgnSpYBoDU5ddJAwNDslTFZT9Zzs58aa\nF3/2Td9rxHQJIXtH+1eQMDy7JEw4TOE30U9v34PdW/69sbWykEwAA8C0YtAym+bQoeBbwvb+\nO3++RU6cUs6OPW3ChO1Gv4ne2YPzB2u65y27REuXHsuDhUfPyJ4lbJv+RHqTeTmecnbs2Kcp\n241+0xz7+xYiyFS8pkuQxf7yo2dkzxI+ugpM+2r6xkQpEu6oXZQm4d9QDx0KxCjpEkJeV8qt\nB/nhg8GzhM2w2re5fEuScJNyJLx3rcGBtq+aRkqXDKp61mKQpZZQZKC9Xm0SllubWamBJLtX\nlxNdg0JseLUgrFZC5dGyVGkkppLwwsRp9HKtqSRcbjcmvFeXIw+hXrNWDpYWwL1sFHNjO3Hl\nEyvb3L9glRe7j8++7FqThIsFXvhu01Q7srQA7mNdQqWWmq53lD+mtN8b1bA6YhihsQgJY7JH\nwtTjhPxzE8++d0gICYtjtWrjFm7vEp5PmSq2QdvjPWUO6YbbdJbpJma9ke/U+oCEnpjHRuvL\nDplywG3nmG5y1uqaTv1wkNATc9nWR3R9pxyfwgLoC4cRKUjoh3m1U2mjUxonzD/d8oCEfliR\nMFbSsSkrgCmBhH6AhMAZSOgJS5tQu0l+6JTjU1gAEwIJPaHckVu5J1SMKbKQMAU+IwsJvcFU\nD/vvF6xPRtsuJIzP7JR7bmuHF/ii0BiKamnAJqEx6gEJ4zNFeWUAaklPy6SO9WQCUmYM5cCE\n8Yb3NCBhOowo2/bE4hW/lvchoV8iSGhuGRJGZ0cIlvS0vQ8J/SLDsnyS9JaE/jc2ZQZwF9sh\nWFpifR8SemZyL9jN8yFhesYoQ8I8mbvnffYo2oTJETFdrO5AwrSY5Z73eil6RzPAGJCaL1ce\nt96HhKEJUWvEOGE2LA9EcKue6B1NgSZhiNF7SJgnGCfMB9GGZ8tnRx8JxKeWALqy/4QLCYPD\nxib84dtR7t9+EqoJoBNHTriQMDijgcHGFCBhjhw54ULCCKhjFRwS1sChEy4kjICQDxJWAyTM\nDT0iaBNWACTMDqWHFL2jdYA2YW7IK30xTlgJ6B3Nj4D3uYCEeYJxwoqAhNSBhOSBhNSBhOSB\nhNSBhOSBhNSBhOSBhNSBhMHBD8KAGfuuCIWEfsBPo4EZe++NAAn9EGaq2jyF+FQSwBAYBwUk\nDEugy5dsScSmjgCGwDwoIGFYICGYAQnjAgnBDEgYGbQJwQy0CeOC3lEwA72jscE4IZiBccKy\ngITUgYTkgYTUgYTkgYTUgYTkgYTUgYTkgYTUgYTkgYTUSSgh8EToSCGAoVncxTHjuYVzZty/\nBa01s8fXV/O0HSKbyeqAgITUgYQum8nqgICE1IGELpvJ6oCAhNSBhC6byeqAgITUgYQum8nq\ngICE1IGELpvJ6oCAhNSBhC6byeqAgITUgYQum8nqgICE1IGELpvJ6oCAhNSBhC6byeqAgITU\ngYQum8nqgICE1IGELpsp+IAAgAaQEIDEQEIAEgMJAUgMJAQgMZAQgMRAQgASAwkBSAwkBCAx\nkBCAxEBCABIDCQFIDCQEIDGQEIDEQEIAEgMJAUgMJAQgMXlJ+Lqw5tEeXOnROKzknpzgn9uu\n+9wZu38d08wb90CYvHwclqeiO9F28fp4yM7KEZOVhI/+B6SaYzvu2q90iZWcoG2cdt37TJp5\n4x4Ik4+P34E7FV1J02/Gh4XLR0xOEn7Yve3OgvcjK/1jzYd/GvYvTnITN7cDpfnltr2xh1Oa\nWeMeCJPfRs4flueiO/HoNvBgt9P5WTticpLwNmTmWAQe7P17/GPPOMkJ/lZ+83F1tU6/ljUu\naeaNeyAMXuzqQcJT0ZU0rPWwlY6VIyYnCUeOfeMb6xpYH/dzldMO/joeKHcvFZssOR0Iwe88\n5e9nif1sycNJc+2IyU/Cll2PfJydPOMdTG7kyr5OKV4YfzZ9Rak4zgZi4uNLHe4aXZMHe53e\nxtoRk5+Er75as5uzsT+Y3MCT/bmlyNit7y5wWDV3vEnoayMdTtE1+NUjzzfhV4+Y7CT8Nsfq\nMydjfzS5nr7O5Shh1zFzP99wyo8cJXSKrsnr1pyO1/oRk5uEbXOw/nAu9oeT67l0Hd+OEnZt\nwq+XjvzMyFBCt+hauJ+tj64fMTlI2I/EjBm5Hj08m1OxP5xcx72v4zhKqP4pinOB0PG0f5yi\na+Nsd/bGEZPD0SAl/F6uRyeTDJ1yX6dOOYfkOtjE4VU9dZznyJlAmHjZP47RtXIyQxtHTFZH\nw9uhL+vZn2TeLm1nl+Q6Tkg45Pbrpc8uM04EYoYPCV2jqzOME55tPxCS0OnYdJ+ocVIFx/HF\nS9t1zPydSTlP/M2Y8SKhpxNdP2OmvXkYo8i8Oiq4OxUwl34Vhz3ultyE23pP19zmj3Mg5niQ\n8GR0J5oI3yonCd1qeW0/eT9acsrqTqu9r265zR/nQMzxIOHJ6Ep+3+rioxwkIiEAVQIJAUgM\nJAQgMZAQgMRAQgASAwkBSAwkBCAxkBCAxEBCABIDCQFIDCQEIDGQEIDEQEIAEgMJAUgMJAQg\nMZAQgMRAQgASAwkBSAwkBCAxkBCAxEBCABIDCQFIDCQEIDGQEIDEQEIAEgMJAUgMJAQgMZAQ\ngMRAQgASAwkBSAwkBCAxkBCAxEBCABJTsYTDD7lex99Yf11Y82jT5ggco5QIVi8hY5/uxaN/\n2j9jas0AAADqSURBVJCMYbWUEsGqJeweH+z6e/yw+y96L3ZPnCdwhFIiWL2Ew58bU98CNCgl\nghTz7AkRwmb2FqBBKRGkmGdPiMrMa3qn7Ss2gAqlRLBqCXse8p0Xe6fLDjhMKRGEhNePeOPb\n3FLmBxyllAhWLWH3+G7YOM7UNhSrMjVTSgSrl5B/RDPiekmYGeBAKRGEhOPf7+X6TZkbcJxS\nIggJ276H+02yW61ySolg9RK216537Us3ghVTSgSrlpBN0w3v4kXqTIEDlBJBinn2xBjAfuI9\nIxzCeiklghTzDEBRQEIAEgMJAUgMJAQgMf8DGd7OPEIAAHcAAAAASUVORK5CYII=",
      "text/plain": [
       "Plot with title \"LAD - t errors\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4QAAAGkCAMAAABeoxsWAAAAM1BMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD/AAD///89ODILAAAACXBIWXMAABJ0\nAAASdAHeZh94AAAgAElEQVR4nO2dDZuiIBRGscxpmnL9/7928xtQAfUiXHjPsztTSUDeOQiI\nJhoAQFBE6AoAkDuQEIDAQEIAAgMJAQgMDwmF0Or5ehRC3Kr3VeWBnazswer7WiVt77hXn1D1\niYeIqyah78JqjOHPNeWBvazswS5gypOeV6D6xEPEVZPQduFzjuDfFeWB3Sz34G8Xr995u+cY\n2uoTERFXTULbhbdvv6Zums9diPKK8sBulnvwG6xvB+aubq9f31gWQeoTERFXTULbhePTeuX1\n1zfYj3Gs+Hp8X3q8pq2fWzsqaZM9b+L2bYGfhbgPDfFvKdph5melvLoqRDEMXuRc+gdbxdQ/\n36qI8rfJkcUf/aeV7WvcR99+M3VI84goSwm/cSvX+jCia2xb+pjdh/5OOW69dU1xOx/Qvf6p\n5u7QmLZ7ppb3KaRuk5xL/2CjmPFdU+OfFQsJq3YA/zNNzczbX9+/dlM2GUSUpYTdmLB4/OoT\na2KiC2s5PS3nrb9SsmLe+s3yXnd/KuWivHHfFyu5/G4W8+g21t94Pn3tlohZSPjdiXXbdSn0\n7d/XbqZsMogoSwmnvXR76cmK7yuPPvm3iRXPb+C/zW/X4RF9VMYHrcm3d/eraaaOUv9MKa+P\nZt2HQM9luxjR52j8G0sXPWKv/vhRjn1PafvCVzWbDCLKU8KmHc93VFqydqcNQ8XH2GJVfTsq\npD+AP+XXoiTltbJtw7sXSy2X7sFWMW1r+7hk+j1GVprNdl+8pC7eVlItmwwiylTCb7f+99F1\nKp5TZ2RONu33uk+qPR+Sqb/aHKu7WAnZ1Dfpei9yLnVjKOanf0+mHmoRm/qhxbzX1pJOodS2\npR3RUxI+2xmSK/7INhrLT9n1DLYlVN4+PV8N2e94bDWEbC2X7WLGFQXFNWtC4kKLmHRmtz/K\nzNuV3p1VQuX1VCJ6TMK+VsMcUmVLfR41pGNrOm7YlnBq0IrGErL2VPLt8XyvhKxY6d00SmHr\nxTT1b7+Hcpwe1SScfBCDcvN2ZXbUKmGSET0hYdWfMa8umP1TQ/qYwvYRq7uz/10uu/ZKMuXX\neLZqJWSlciJrEbLNYjpeD+3PMRPUT/0nZNRx2912nnD+nWpET0hY9O1FfcHsn7Ybvk8fbY/g\nVajnmNSQrUxyKcmUX8OTtXbz26QWf92v+0oum8XcphHGFStCYkONWDWvV3v2Xadx+19p3D95\nRPSEhHPd6aqzWd5MI5/HmVZgjMnk39PZ2t5UY8ju3d9Ga/UiZNNZJbUJnx5sFNNOyn8a6fx0\nVqgRE/MAop728IRhViGPiJ6Q8DHW3X9Lr0k476XipSVTft/lXWkO2dhhKrrAqCF7DduqtVy2\nixmH8TkOCdWI/codllI9vW6+iiKPiB6VsPx5vrouRl1d0NLrEvbXE34rUevJlN99MmkJoPJA\n/fV+tItw3p/u3JEasu9nnKeBV0K2UUw/erjnuF5Gi5gy7HuNS8T6v+eqNuSSSUSPSjgJ8f2g\nxt0IADBzcDz3fj+fZdlNz5ibMgCABe+TKgIQ4TtSCKBvNnex9xj6LiAXgkkYqNz0gITsgYTc\n8SmhuaOEGBIBCblzrYROHWGwD0jIHXRH2QMJuQMJ2QMJuROFhOibniEGCRHBM/iR0OWi3qkA\nZdE32E14CRHBcxBLqCyjNi8dFeoDhPAoEUgYtB788SGh20W9YuM32ElwCRHBk/iQ0O2iXkhI\nBCTkjg8J3S7qhYREEO+4umovA/1pby1tvq07JKTCh4RuF/ViTEgE7Z5rb+ne1MMV5saLVDEm\npIJcQueLejE7SgTtnnuIsv7+aG/V83kYI4jZUSrIJZyWpNku6sV5QiJo91130xcxjupNfRmc\nJ6SC+jyh+0W9iBoR1BI27dSa9ETdisW/9ESxYgacgbo7+m7vG9Z94VhtHBQigFRwkhCN7yq0\nO+UtiurdlMW7+3Id4w3PSMvNGT4SYvS/AfEueRVzh/PnwnIzhpGEh96VAeS75PfRfUNE+WP+\nehOEgooYJVztduKM8BbBV8yAk8Qn4Ua3ExJuAQm5E6GE61sh4RaQkDvRSbgpG8aEG0BC7jCS\nELOj60BC7vCREOcJN4CE3IlOQnQ79wIJuROhhOh27gMScic+CdHt3Akk5E6MEoJdQELuQEL2\nQELuxCkhOqQ7gITciVFC16kZqNoBCbkTpYROFcAs6gAk5E6EEjquEsX5xAFIyB22EmJB9wg/\nCTGOUCGX8O+n7C7KLqu/kyVDQje4SYhxhA6xhPVNuhuX8daxZ8eEkHCEnYSn3p0ixBJWovjt\nbtTVfF6F481/l1ucmkrRJ0UwmUmI1nMBsYRFf6+8jrfjbfDXtjmdoOjugomODSTkDrGEihFu\nXwhznEHV3KPJQUKpVYWEC+I8ErqBcHbEL6E6vsCYUId+TPjqb5R3YkzoCiTsYCChkh6zozrU\npyju0uzozfULYQ4CCTuil3ARJ0ynqdCfJ6y684RF+XP0PKE76Ni08JMQqES4YmZHFujYNJCQ\nP6wlRMemJXoJ0WOxwFxCwEJC9FiM+JTQ93lC0BGzhGNPBT0WE9dKeOaLXrfekX1845UQR0A3\nuHRHt+KJOMcsoWO63GEj4cZ7EOd4JcSsqCNMJNyKJ+IMCfnjR8LnTYjS9H3newuAhNtAQu4Q\nS9iPzoa1a8aloxdJmMG0TbQSYqzgiA8JK1HVTfOpxPNQyYbUe8eEWUzbRCyhmH8kHoRT+JCw\nEN3K7VrcDpW8nfHu2dEsmuJ4JRzly6ItPIEPCcfdTXuy3naecGV7HoOSmCWUU6YdhTP4kPAx\nSujnot5VG1cbW0gYRbl5hOEE5BKWP8+X+P0+rCsvF/VudG1WG9s8og8JuUMu4bQkTYjCx0W9\n612bjThn0Q+ChNyhPk/4fj+fZdlNz1RGB49avv7urZdzmBGIXsI82sIT+DlZf6pkp3e5djwz\nmBtnIGEObeEJEpEw58Y2fgmzaAtPwE3CLdsybmw5SAhM8JNwy7ZsG1tIyB12Eqq2ZWueBCTk\nDkMJx/cLkXEfVAIScoerhJJ/uf8xQELusJWw+z/N0uTcLYWE3GEq4XgI7H/n3S2FhNxJQ8L+\ntUw1hITc4Szh0B8Vo4kizwtI+UiYW2RcYSqhdAicvyq0+ze+lk+4uUiY96DBBFsJR9fEOEnT\nzIfF/tVcws1GwkPvygGuEip9m/Eo2Mw/RD7hZiLh1qrfOUE2zaYOXwnlrKbxoHpUpCshZlhL\nOJmXc2c1CQnnWLY6ClXC5BtYxhJK44acO6uJSNjIEzJCSEfF9BtYJhIuNRsVlHstSUdqi2Qk\nnGI6/mhyaWDZSKi3h323RR3AJx2pLRKScMx0nDPNpYHlIqE+MphOJ0HC3Ru8l0yR6TRG9FlY\nJPCRUH+/kH6vd1bzwI+E5F8I44ouHSSMt9xhNDiveJI7q1NXNQeIJfT1hTDuFdDyxpgw3nL7\n/qh0kmJxHMzEQh8S+vhCmF0VUC69l58n2bbylXBQcCUjkUPrOeFDQh9fCLOnCmLjeaKNK1sJ\n57nRlS3q77TxIaGfL4Q5T6KNK2cJm42/Ekho3mDOrn2b5y+EOcg89ri4YN94+TwOU2sU5W6O\nD+SJ0+Qhl9D3F8IcRFtemhK0n8d9as3rflSuUkv9eEguoe8vhDmCui4jtYh6kNBpas3vfuxi\nNi+6SC1oCtTnCX1/IcwB5tVss4gp4UFCp6k13ztSyB2X5KIm4+dk/amSfRTUL42au6Td64k0\nrx4k3JwXETKk5W7WR/udIulLODan0iCjeyWZXo4HCZ2m1iAhFZlI2KgXzEgv84daQtepNUhI\nhU8J4zhPONim3XYmodhSS+g6tXbRrlNay1SGECrXSnj9kKKRLBSLFyHhEueptasknMcN6Qwh\nVBLujk5X2a+FDhLyKVdbZpFAyDSSlXBWb/2wm05A05ZQvanehQVfSLoSdj+2z86n07VJWUL1\nhMmFBV+KNwmtf99+d+XGiQklSSKLorhLaAqAOitDXHA0pC+h+RrtBA6IvCU0BkDTLp0hhAqx\nhEI4T4BeJ6F2dmKlFqzjylxCU166hKOw7HsvKsQS/hWRSDiNCc3XTqTQw2EtoTkAi63SffTm\nd7BXkro7Wpfi/uneH7Y7Ot9Obxr6mSrBOYoJS6gfJ6Ubys7vTHhAcfhT/YpuzVNoCadWU+7K\nbFaCcwyTllCZHVVn2dIZKHqYmPncRVlHIKFSztY12imHkEW5tgAszhPKEoqFkizxMjv6I4pX\nPBIq7efyNlBNk2pnhkW57gHQdVOXQnEOoRcJm/fNPli+bK9N12ivRzzdYT2Tcp0CIB/zlL4N\nJNzmEZeE1vEDYxW5S+hSlDo01Po2F1eGHl8SniiZuBgpflvjh/U13lxIVcLlaLBRDNSma7iG\nr8lBwsVPtei+s9pIUeVGmhKuNJ7DuGKxlHTrSpk5q8gjm7qEinfavJpopOUXouFqYaISLn5O\nA0OxTGKYY2XQy8lGQqkHI+Y2tf927TGOUUdqkyQl3Go85d/S1WqGKhn8jIVMJJwGDv3j8Y4z\nY1d0CGTUkdokfQlXzhAODyyDffOmaEhdQr1bMxwQx6iK4cVmfpEdOUiozIeqJU9dmfUqQcJD\nJRMXIw8Jpl+ShIOFIvJxwzZJSqj3IpVVo/K86azn+vwLJDxUMnlBK8sNFQnHoSFTEpVwbT5F\nGgUqzeo0vN88FEYd3QwkbJbN4diuShbyJU0JN45sQh0GzoFVJ06VdzRLm+MiDwnVkClXPMpD\nCp6kKqGlWKUHalzHHX0bm4uEShdmvuZYOSUce6w2yFZCaSTP+2KKTCRU5tYG8dRRvtvtAGIk\nXwm1nwErdY5sJNQLlhfhT+tm5qMjIzKTcBE/ec0hs9AN5Cfh2up7aYomYM0Okp2ESk+mEUJt\nTRmSoYRyo6lPlzYMLcxNQr2tVGPIkSwlnBtNaTjYQEIe5WoVCF6N8+QnYa/afLZJ6Ydy7I/m\nK6GtB8rl2JijhOOYYl5Guj0mZBDHbCW0RIfPKDFDCfvYaddVNPMJCilu6txppGQsoRE+vVVy\nCf9+yu6Puaz+DpZ8DdK0qJAeqbZNpxKbiNtUSLgKozOHxBLWN2lJ2P1YyZcwT6lpL6wm2kwR\nA5BwlXwlrETx++4efV6FqA6VfAlziAyGKWcUQ9d4E0io4HChb2wQS1iI9/T4LYpDJV+CJOH2\nYjVIGGO5RvTTT3HWUoNYQnVEZcwjkjGhScLpNKL8luiAhBKSeXGP5GVSPxJuTmvOE59SbfTU\nQj5QxhpOSDijjeJZKOhjTPjqvhktjjGhsTFcLH7qn//rkC57mn74r+8hIOHM9J2woSuyC+pT\nFHdpdvRWHyqZEIdhgSxh+/vfv9FCaWARdZuat4TaOaXxlG+w+hyB/jxh1Z0nLMqf8OcJXSbI\ntGnSfxON0qRGHNWcJdS7OtK6J0Y4SvgsxO15Ucn0RTgcCgff/v1TLExIQvoIxrBLtK7OfMuS\nGCrnjE3CdymKZ/PjcPKdrGT6IswSSjM0/1QJpz5qFH9wW1jq5i2CEewToQ4aphWG89onScd4\nzbRI+O5iV4lH3XxKQdqSRjImbOboCFXCf2KerYk1fC3muvmLYPh9os1wN9OlvkI2MP6Lfi0S\nPtoZzqo/11CL2868nWZEfLJvzwuhSjidb4ozdCPm2p2K4Ilyr2A6DGpHOyEa+QipP4wOi4TD\nZymlJzvyXqQXMvsyO4i9HDmFJmGcIdMx1/JUBE+UewFiaiXnl6bZ0U5GdcA/PIwPJwl/+16M\n+eQ7WcmXoh4rte5oHFW04SKhjwgG3zurgwX5crREJHy0Y4me+mE++U5W8qWovRRtSMgDW3fU\nVwSD7x95+Ce9Ko0G05CwLuZrfWgPhHHsDjU26sxMFBV0wFxPfxEMv382xnnS9dljAv3Siriw\nSPgd04+BK/a1ovax2K7sPKFLKE3NSBWMe2rGVrejETxbrn/UzqiQHk1XaI9DxKiXHlolPJwx\nTwn1ZWtN3JPbLaFqFsMeWY3S+HN2T0sbG8QSCpVDJV+K2kuZLRTTNHfUk9stOUs4szgbMUrH\n4Bb5jhL+3UVRGddjD+kKbhLqPZrx4qXu6fSjiaW+K7hVzDWC1OUSlWU7islRmg+EWnR91e4s\nNgnf93bRU7/sonCIYV2Ke3ctE4/uaKPHd1p82D9rFhLG16mx1GdvBKnKpSzJPiBQJBxmRoX0\nVE4THRYJ//pFT/fi3dR3twnuXyF+G0YSyozT22KKpCZhjONDc20ORJCkXEocBgR6lKQ1M8bb\nCEWBRcJ7v+hJvJp20ZPbDPfnLsqaqYTDT0lAsTIyjKvq5tociSBFuR4KMlso5hZTGg0t1v/G\n15OxSihPN7kfAX5E8eIo4ezgGD9tcjvKjo25MscieL5cDwUZCpTOSIzPmzl4UlMaYU/Gk4TN\n++awZtM1s+uQJBzip7WekNC5XA8FmSRsFoP7ZvFSM0c4LvxI2K6WYi3h2JYqGyDhjnIpsaqz\nsiBm9ZgXZfz8SXii5IDoA79F0GJsSTOQ0NyJVCfSppe0zmnT8JXQ29VHke2IDnnx77qEEY4p\nbBL6iuCVe8FYc2UirU/dSFGS19JM74gKSKgiJA8baS9oCWIiBwlNzAaqjq09i7EnY5UwRMkx\nsbb0NzYNQ9Umlr0wHug2Op3Ksxh7MpDQij7mjy+MkLD7OcdkQ0Ltts4RAQmNyKee5F8xVT53\nCRchWe2HzjcQurRuTkBCA0Lu5ejtajT1h4R652S1Hxrx1dqQ0IAySaO1q8pRMii+KuF64UIE\naJFYzo6qN9KLI24zkHBkqZSy2lDrjYrFfE0wIOGSxaUxUd9GDxL2rM23COWwN77YDJPh07/Q\n0FbB/ZRGBB/dGf2OsnHVHRL2rM23TFPe6nnCZlpVOvwP3C2lLdz9suy4AmgBEoYpYA/aFOj4\nbD63NDyfB4nT/+DnLIiLNl6W7W3phm9WJIznE0DCji0J5b82dfnT2BVdLh2+GvKiHS/LjiqA\nNrTbOqvRDG0jJOwwSqhsm8zTDoHh2lb6Mt0uy44qgDZUB4UUzeA9GUg4sno8U14UUj90GifK\nioaKpo8SXS7LjiuANpTeqNzqBu/JQMKR1dlR5feYQkhPRAzR9FKiw2XZcQXQiuKgtOhifBCm\nVuayM5NwrS+50g8dXZO7pvLrTYDP5alA62XZsQXQRBtcZdnaqoShBoeQcBtFK7EUb4xcMx8h\nm3QkjLbc/QwhUpaOLtvOcINDSGhAOwQqs6PypvULgK8CEtpYGyhM0ZxnaFZSXYMfCZ/fEUX5\nOlhyPMht49hdWZFQGR0mMiaMuNzdbDSPQjoNHLInQy5h/7nu/dS++U6zLGK4uEpNrD2RlnMn\nMjsac7m7WdNrebOExCSsRPulB59KPA+VHClz/0WoXdNhaakINbCHhBaWeokhhoqKi1SX4UPC\nQnRfeVCL26GSo2VeyS3HT15DE6ZamZW7H32gMJ5gUrsuyYwJlc9lPi7wiaHCIqCLebarya3c\n/egDBXk+bf4YsprX4kPCxyih8ZsP+MRQZtFpgYQcUNRSrhNd9FIDjOvJJSx/nq9u/W9dmWdm\nOMVwRpJwGs/Pk2xBq5RLuafZlDBUl5RcwmnNs+3b8HjGUDkpMZxnmpdzB61SLuWeRqjNpjTH\nHWhyhljC5v1+Psuym56xfC8s0xiqnU/RhFYQEu5H6XWOh8X5FfYSEpQcN4uLJwIrCAn3I6Tv\nRFdVHBNISS/4mJBwN9rFE8ttV9fn8hLDlruXtZhMq2XULqk+JrxongYSHmNVwsjm1hItdx9u\nMZmH+krqi+ZpfEqY5HnCkbX4RDa3lmi5+3CLydymhlhEc62EbO8TtERvM8X83aLSdNslNbmk\nlHjK3YXVI+WUvT64sL2ZCHRH9yFU76aHzTjcH828rmsKCQ1YPJIXBAs1XMrdhfw2qJBwD1tm\njdPcvYLbTaufSl1QRkzl7mJbQuXys2G+VN7aTJEWS0P9VNJ9g/eSI2bDrKn/Ofzcnj31V6nr\n4RHA1bH7YNV8qFukEvPr+i29fNVxzwYnnoW4GS9k4hJDhS2zpPYUEkbGsvOimNU0sl/atOjY\nSx1e8vaBqSV8l6J4Nj9dW3M/VnK8WCUc5mUgYUzox7DOvmYauM8rnhZXcA9plY6Olxru3mDk\n3V9SLx518ynTuqi3MZgl9UOn5nU9pcdaXQ3DALbosRrH8cryQ9YSPtorJ6r+GqbkLupdmjUt\nxJ9mR6fXMDsaKZqEkoOyhtJP/3eUJZZw+AMspSe0BYRFW+arzqDpS7lxnjBK5BHeNDuj3DSh\naZax9dqgepHwt++HJnlR72JBhdjYehmQcB9yR0W6eknIdjZyLBfnL7zUaNcGI492NNhTP1K8\nqFfiwrkXI5DQBaF3YKa+y/iiYRLUc+NKLGFdSJ/VeCBkFsMVDBJeekSEhHb0cYQ0GSqP5xt5\nYvvC2u3eYKEa1SvMtx1lFcNVtudKL5uTWS8/8XKPYJxaUYeFMQ0ovNeFUwzX2YrsdWcnri0o\nknIP4DZykFb+XgskPM7GEe/qsSIktOIakkDX9kDCMywWY8wjfUgYEbHMoW0ACenQ1j5Bwni4\neISwE0hIhzrHjTFhRLjNlYW60hwSkqFPgF9d7tXwCqBdMC1sFxoJCcmQuqGOAaSJMySkQRlM\nSGf0ryl53wbvJbNi4/KXxsEwqgMmJCRBKBFRzuNfUPTODd5LZoS2EEP+aYmfmNcqnq7E6Rx4\nlesJoZrXz3OL5orPCQnPoHqkrLxXtizeJ3VZT+8HSEiCvHh0iuR8VPRb9N4N3kvmw+JkxMYd\nnVffSDaLCglpkBaPqutIIWHMbKtmlpA2xJCQhvnwJ9RZGtMHJZm5gYQnOCUh2YgDElIh5vGg\nPkez9YYp1alyd2+gIoUYbreUxjZ0nLsRFBGEhGRMVzWN1/HaHPM9oICEDmxHyRy/yUKKnQAJ\n6RjFGzsp5ghZRv7upe7eQEUaMdyO0nLL4uJuoirQZMOmXL+MpylcbrkNCbmxdnE3TcZUGTEp\n1y/j9LZLfCAhN4jGD1sZX0+aAdzVSmJMyAuqVnM756tJNIBjf9RlvBDr7OjfT9nNK5XV38GS\n0wQSMkFZtuaQOr7zhPVN+hbQ5L6L4gyQkA3zuQmhv+6pwN0bjFSi+H13jz6vIvH7ju4EY0I+\nyKspptcoZ7TXituzwUgh3tPjd4p34D6OtyBCQmqUiwqnF6Wf1OXt3mDObtF27C85WTx1ZyAh\nNWtr6/0NJ3AkTABISI1Y6bYwkvA7Jnx9ukcYE14FJCRHLNcUMpKwuUuzo7falDLhGF4LJCRn\nbfzOZ0zYNH9Vd56wKH9wnvAaIKEH1lb+NgsxqQrbvcF7yWAfkPAitibWTk+4QUL2QMKgEBwh\nIWEYcBVFKhCMFX1KiPOEW+B6wmSgmDW9VkIhc74AtqxMgZ/KLAgZB3D7ls/Hstu9gYqcY0h6\nKISEF6MEb+iNnoolJAwA7Rc3QcKLUYM39mpOaAgJA7C6QPh4bhSZMCo3NGoHVAwGnjkWepLw\nWYjb82DJGTC1nkSZBSHXAMoSSt8LeyKa1BK+S1E8mx9c1GtkuOUojoQcUSRs5vHg8WgSS/ju\n7KvEo24+pTAeC3ONYYfz7RNc8iLJRcfel8k2gPOYsHs03B/xRJNKLOGjvXKi6q9hqsXtUMkZ\nIITchJ7NjCKTGee+TLYBnGdHewmncB7OcPcGc3Z9tUrpCW0BaSDGHilRZoS492UyDqD6VTGn\nm1QvEv72scNFvVtQnOFdZEaDe18m5wCOTP3RU20qeXf0MV5EWD9wUe8m8Z4nNPdlsORJZbo/\n4tAxPZjL7g1G6kK60bvxQJi5hNGumHHvy2QdwAn5TFMkEn57MmPcCuNxMPsY0h1JqLujrn2Z\nzAMoM5+3P/H2PRuoQAyJoN2R7n0ZBHBCmi89slcgIXuod6RrXwYBlBDqhOnON+/eQAViSARW\nzETBiYEhJAwLwdAQEkYBJGQKySQpJIwCSMgUktOFkDAO1oPp0tWBhAGhWTgDCeNg9Y7BTl0d\nSBgQSJgWK3cMln4a3rd7AxWIISRMnXFht2VXQcKQYEyYNo539IKEIcHsaNo43tELEoYF5wmT\nZnTQvK8gIXsgYbyMXR1ImDiQMF6mW19YUu3dQAViSAQkjBiMCfMAEkYMZkfzABJGjcPUGyRk\nDySMHZuHkJA9kDBuph7ppouQkD2QMG6GuRmDi34kfN6EKF8udQOngYRRoy4QXr1ZN7GEfe73\n/r6UuEXJJUDCqBl2k/alomtJNt+7t8T2bZWo6qb5VPhCmEuAhFGjSqi62Kw8cdpgLrF9WyG6\nO1fiC2GuARJGjXp37sskdPsiWsSQCEgYMeP3VOh3RPQu4WOUEPeOvQJIGDH9UbCfFG2ahYvN\nyjOXDeYyRfnzfInf78O6whfCXAIkjBflwCe7uJpq8+07C52/sEeIojYmPVQAWAAJ42Wt9+n/\nPOH7/XyWZTc9UxkdRAypgITxsjoE3EzlvoEKxJAISBgxa0PAjUR7NlCBGBIBCSMGlzLlASSM\nmsCXMuE84SVAQu5cKyG+8twDkJA76I6yBxJyBxKyBxJyBxKyBxJyh1rC+iHEfbicFxMzlwAJ\nuUMsYV10cy5lnwUk9IQyrQUJGeIWwWO7uLuQt34W974gYzUOFQAWJ4AhITtcI3hsFxf92z7F\n7QMJfaEthYKE7HCN4InrCb/U9zsk9IS+KBgScsM5gsd28U2Ml07c7pDQD5CQO54lfIrH8Ogj\n7pDQC5CQO54lbKrJvJdlZRpieBSMCbnjd0zYNO9yfPR5QEIvYHaUO35nR/dUxHcBCYPzhNzx\neZ5wTzV8F5ALkJA7kJA9kJA7kJA9kJA7kJA9kJA7kJA9kJA7kJA9kJA7ASUERPiOFALom81d\nfGU8dxfHObnv7K/HSw0zz9RzzhTFcU4OCZFp8JwpiuOcHBIi0+A5UxTHOTkkRKbBc6YojnNy\nSDQGmCQAAAVtSURBVIhMg+dMURzn5JAQmQbPmaI4zskhITINnjNFcZyTQ0JkGjxniuI4J4eE\nyDR4zhTFcU4OCZFp8JwpiuOcHBIi0+A5UxTHOTkkRKbBcwYAOAEJAQgMJAQgMJAQgMBAQgAC\nAwkBCAwkBCAwkBCAwEBCAAIDCQEIDCQEIDCQEIDAQEIAAgMJAQgMJAQgMJAQgMBAQgACc4WE\nVSGKqpbKNH9HTdM81W36+83JLbk/b1pm5tz15Jbc64cQj7dz5RfJrbsmGLYgHORJ/1kXESZg\nGVdKLoj3vfu7uk3P37a/tLe6TX+/Obkl96rbWMwxMueuJ7fVvei2StEyZ68nt+6aYNiCcJA3\n/WddRJiCRVxJ8R/vP1G8m3ch/sYX3qI0vuGbVq7V4v3m5Obc3+JRt83vwy33RXJL3as2ZSWl\nMWe/SG7bNcGwBeEgWuxIstRDRsEiULT4l7ASr+/PX/EzvvCcH67xFHclMov3m5Obcy/7pPM7\nzLkvklvqXohaSW7JfpHckn04LEE4iB47ChYho2ARKFr8S1iKT6O08U/xNNaoUj/t4v3m5Jbc\nhzdN77Dkrid3y72YHrpkLyV3yj4ETrtpN3rsKLP2kK8UKOKMPeUrlaA3TaV4Pb5D5630b20P\nLt5vTm7JvaMWd8fc9eQuuVeSSQ7Zy8ldsg+Cy27ajx47OqSQkVF5ayGDSNhh2Et7JNQ32XNv\nDzgv99yV5Pbcf4WQJLJmryZ3qXwQ/EjoJccOKWREqIGiJYCEQvx+mypTu3JGQnvuzacopeS2\n3LXkttyfZSGNnKzZ68mtlQ8DMwnlkBGhBoqWABL21Ib57jMS2nOvC+lIY89dSW7P/ctjlsjl\nj/ehO2fJPgS8JFwJGQWLQFHhX8JiI36uVm2935yTIfld/hO3535fMcL8p1PPI3h79kpyp+xD\n4PI5DuHlo66FjIBloIi4anb0s5hYcx/krb/fnNNm7p/b/bMjdy25LXd9s73ya7nFJ6HT5ziC\nh4+6HjIKvE3l+slW4qcbI7/mYW1/0sXZqsX7zcktub+0WQ9L7npyS+7j5qklNme/SG7fNYGw\nB+Eg9H/XesgoWASKlgArZqo2lHVlmL/atWJGS27O/aNHyJz7Irml7t3Kirqcxw4OK2bk5PZd\nEwhPK2Y8SLgIGQWLQNFyQcfnNk+7d7u87hfiGRrVMTL975vrGQ2H3B9iXp3pkPsiua3uhfZZ\nLZXXk9t3TSisQTgIuYRyyOgoPH38ngskrLsF+H1pYnzhZmpUVAml97sm38xdSCFyyH09uanu\n82anyq8lN+6aUFiDcBByCYUXCRuvcYluCgCA3ICEAAQGEgIQGEgIQGAgIQCBgYQABAYSAhAY\nSAhAYCAhAIGBhAAEBhICEBhICEBgICEAgYGEAAQGEgIQGEgIQGAgIQCBgYQABAYSAhAYSAhA\nYCAhAIGBhAAEBhICEBhICEBgICEAgYGEAAQGEgIQGEgIQGAgIQCBgYQABAYSAhAYSAhAYCAh\nAIHJWML+C13vw9ewP2+iqOqwNQL7SCWC2UsoxLt9UnUPC5YxzJZUIpi1hO3PSty/P9/i8Y3e\nUzwC1wnsIZUIZi9h/6sU8kuAB6lEkGOdiRhDWCxeAjxIJYIc60zE2Jl5Tq/UXccGcCGVCGYt\nYUc1v/IUr3DVAbtJJYKQ8P4eX/gUZcj6gL2kEsGsJWx/vgoxnGeqC45dmZxJJYLZS9i8x2HE\n/RawMuAAqUQQEg6/P7f7J2RtwH5SiSAkrLsZ7hfLabXMSSWC2UtY39vZtQ/fCGZMKhHMWkIx\nLTd8jE9CVwrsIJUIcqwzEUMAu4X3gnEI8yWVCHKsMwBJAQkBCAwkBCAwkBCAwPwHtxlMj8Au\n63IAAAAASUVORK5CYII=",
      "text/plain": [
       "Plot with title \"LAD - normal errors\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "par(mar=c(4,4,1,1),mfcol=c(1,2))\n",
    "plot(ls_yt[,2],ls_yt[,3],xlab=\"B2\",ylab=\"B3\",main=\"LS - t errors\")\n",
    "points(2, 4,type='p',col='red', cex= 1.5, pch=19)\n",
    "plot(lad_yt[,2],lad_yt[,3],xlab=\"B2\",ylab=\"B3\",main=\"LAD - t errors\")\n",
    "points(2, 4,type='p',col='red', cex= 1.5, pch=19)\n",
    "plot(ls_yn[,2],ls_yn[,3],xlab=\"B2\",ylab=\"B3\",main=\"LS-normal errors\")\n",
    "points(2, 4,type='p',col='red', cex= 1.5, pch=19)\n",
    "plot(lad_yn[,2],lad_yn[,3],xlab=\"B2\",ylab=\"B3\",main=\"LAD - normal errors\")\n",
    "points(2, 4,type='p',col='red', cex= 1.5, pch=19)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(requires internet connection to view)\n",
    "\n",
    "<img src=\"https://longjp.github.io/statcomp/hw/betaT.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(requires internet connection to view)\n",
    "\n",
    "<img src=\"https://longjp.github.io/statcomp/hw/betaNorm.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 2\n",
    "\n",
    "Compare the speed of $\\beta_{LS}$ to $\\beta_{LAD}$ using a package such as benchmark in R or some of the tools [here](https://jakevdp.github.io/PythonDataScienceHandbook/01.07-timing-and-profiling.html) in python. Empirically determine the computational complexity (in $n$) of $\\beta_{LAD}$ and $\\beta_{LS}$ (we know from theoretical analysis that $\\widehat{\\beta}_{LS}$ should be linear in $n$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Installing package into 'C:/Users/Kanchan/Documents/R/win-library/3.4'\n",
      "(as 'lib' is unspecified)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "package 'rbenchmark' successfully unpacked and MD5 sums checked\n",
      "\n",
      "The downloaded binary packages are in\n",
      "\tC:\\Users\\Kanchan\\AppData\\Local\\Temp\\Rtmp04U1sr\\downloaded_packages\n"
     ]
    }
   ],
   "source": [
    "install.packages('rbenchmark', repos = 'https://cran.r-project.org')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>test</th><th scope=col>elapsed</th><th scope=col>relative</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>4</th><td>Least absolute deviation - n dist</td><td> 6.39                            </td><td>1.051                            </td></tr>\n",
       "\t<tr><th scope=row>3</th><td>Least absolute deviation - t dist</td><td> 6.08                            </td><td>1.000                            </td></tr>\n",
       "\t<tr><th scope=row>2</th><td>Least square estimation - n dist </td><td>10.27                            </td><td>1.689                            </td></tr>\n",
       "\t<tr><th scope=row>1</th><td>Least square estimation - t dist </td><td>11.61                            </td><td>1.910                            </td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lll}\n",
       "  & test & elapsed & relative\\\\\n",
       "\\hline\n",
       "\t4 & Least absolute deviation - n dist &  6.39                             & 1.051                            \\\\\n",
       "\t3 & Least absolute deviation - t dist &  6.08                             & 1.000                            \\\\\n",
       "\t2 & Least square estimation - n dist  & 10.27                             & 1.689                            \\\\\n",
       "\t1 & Least square estimation - t dist  & 11.61                             & 1.910                            \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| <!--/--> | test | elapsed | relative | \n",
       "|---|---|---|---|\n",
       "| 4 | Least absolute deviation - n dist |  6.39                             | 1.051                             | \n",
       "| 3 | Least absolute deviation - t dist |  6.08                             | 1.000                             | \n",
       "| 2 | Least square estimation - n dist  | 10.27                             | 1.689                             | \n",
       "| 1 | Least square estimation - t dist  | 11.61                             | 1.910                             | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "  test                              elapsed relative\n",
       "4 Least absolute deviation - n dist  6.39   1.051   \n",
       "3 Least absolute deviation - t dist  6.08   1.000   \n",
       "2 Least square estimation - n dist  10.27   1.689   \n",
       "1 Least square estimation - t dist  11.61   1.910   "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "library(rbenchmark)\n",
    "benchmark( \"Least square estimation - t dist\" = { for (i in 1:N)\n",
    "                            {Y_t[, i] = X%*%beta + rt(n,df=2)\n",
    "                             ls_yt[i, ] = (inv(M))%*%Xtrans%*%Y_t[, i] }\n",
    "                       },\n",
    "           \"Least square estimation - n dist\" = { for (i in 1:N)\n",
    "                            {Y_n[, i] = X%*%beta + rnorm(n)\n",
    "                            ls_yn[i, ] = (inv(M))%*%Xtrans%*%Y_n[, i] }\n",
    "                       },\n",
    "           \"Least absolute deviation - t dist\" = { for (i in 1:N)\n",
    "                             {lad1 = l1fit(X1, Y_t[, i])\n",
    "                              lad_yt[i, ] = lad1$coefficients}\n",
    "                       },\n",
    "           \"Least absolute deviation - n dist\" = { for (i in 1:N)\n",
    "                             {lad2 = l1fit(X1, Y_n[, i])\n",
    "                              lad_yn[i, ] = lad2$coefficients}\n",
    "                       },\n",
    "          columns = c(\"test\", \"elapsed\", \"relative\")\n",
    "    \n",
    "  )  "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}