LinearRegressionGeneralized.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from sklearn.linear_model import LinearRegression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True theta : [0.88971132 0.05193854 0.6667153 ]\n",
      "X.shape: (1000, 2)\n",
      "theta.shape: (3,)\n",
      "yTrue.shape: (1000,)\n",
      "yNoisy.shape:(1000,)\n",
      "(3, 1000)\n"
     ]
    }
   ],
   "source": [
    "m = 2\n",
    "Npoints = 1000\n",
    "X = np.random.randn(Npoints*m).reshape(Npoints,m)\n",
    "theta = np.random.rand(m+1)\n",
    "print(f\"True theta : {theta}\")\n",
    "yTrue = X.dot(theta[:m])+theta[m]\n",
    "noiseAmp = 2\n",
    "yNoisy = yTrue+noiseAmp*np.random.randn(Npoints)\n",
    "print(f\"X.shape: {X.shape}\")\n",
    "print(f\"theta.shape: {theta.shape}\")\n",
    "print(f\"yTrue.shape: {yTrue.shape}\")\n",
    "print(f\"yNoisy.shape:{yNoisy.shape}\")\n",
    "\n",
    "Xtest = np.concatenate((np.ones((Npoints,1)),X),axis=1)\n",
    "Xtest = Xtest.T\n",
    "print(Xtest.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def computeGradientGeneralized(X,y,y_estimate):\n",
    "    '''Computes the gradient for the Quadratic loss function. X is a (1+m)xn augmented numPy array with n observations and m features\n",
    "       X = [1 | X_input].T\n",
    "       delta_y = y-y_estimate , is an nx1 numPy array of the error in estimate at the current iteration.\n",
    "       returns the (1+m)x1 gradient at the current iteration\n",
    "    '''\n",
    "    n = X.shape[1]\n",
    "    delta_y = (y-y_estimate)    \n",
    "    return -(2/n)*X.dot(delta_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def GradientDescentGeneralized(X,y,start,step_size,tolerance,Niter):    \n",
    "    ''' Will perform gradient descent on X (nxm), using target y (nx1) \n",
    "    '''\n",
    "    X = np.concatenate((np.ones((Npoints,1)),X),axis=1)\n",
    "    X = X.T    \n",
    "    print(f\"X.shape = {X.shape}\")\n",
    "    gradNorm = 10000\n",
    "    i = 0\n",
    "    current_theta_estimate = start\n",
    "    norms = []\n",
    "    n = X.shape[0]            \n",
    "    while i<Niter and gradNorm > tolerance:                \n",
    "        y_estimate = current_theta_estimate.dot(X)\n",
    "        grad = computeGradientGeneralized(X,y,y_estimate)\n",
    "        gradNorm = np.sqrt((grad.T).dot(grad))\n",
    "        norms.append(gradNorm)        \n",
    "        current_theta_estimate = current_theta_estimate - step_size*grad\n",
    "        i+=1\n",
    "    plt.plot(np.array(list(range(1,i+1))),norms)\n",
    "    plt.xlabel('Iterations')\n",
    "    plt.ylabel('||gradient||')\n",
    "    plt.grid(True)\n",
    "    plt.title(\"Convergence of gradient descent\")\n",
    "    print(f\"Started at location {start}\")\n",
    "    print(f\"Finished with {i} iterations\")\n",
    "    print(f\"Finished with gradient norm = {gradNorm}\")\n",
    "    print(f\"Finished at location {current_theta_estimate}\")\n",
    "    return current_theta_estimate\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def VanillaLinearRegressionGeneralized(X,y):    \n",
    "    '''Assumes that the input matrix is an nx1 numPy array with n observations and y is an nx1 target vector \n",
    "    '''\n",
    "    n,m = X.shape\n",
    "    print(f\"(n,m)=({n},{m})\")\n",
    "    theta_eatimate_start = np.random.rand(m+1)\n",
    "    print(f\"theta_estimate_start: {theta_eatimate_start}\")\n",
    "    learning_rate = 1e-1\n",
    "    err_tolerance = 1e-6\n",
    "    MaxIter = 50000\n",
    "    estimated_coeff = GradientDescentGeneralized(X,y,theta_eatimate_start,learning_rate,err_tolerance,MaxIter)\n",
    "    return estimated_coeff\n",
    "\n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "\n",
    "    \n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(n,m)=(1000,2)\n",
      "theta_estimate_start: [0.65370528 0.70731529 0.36034643]\n",
      "X.shape = (3, 1000)\n",
      "Started at location [0.65370528 0.70731529 0.36034643]\n",
      "Finished with 65 iterations\n",
      "Finished with gradient norm = 9.807572611845027e-07\n",
      "Finished at location [0.61444751 0.76669669 0.01412282]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZwU9Z3/8den5z5guAeG4TQoggcKiPEKJJpgEjWHrpqsG3djWPcX1ty7uNmYxCSbO5ts4iYx5tw1IcZEw7okmlVGY+IBeICAGECQEeUGGY5hjs/vj6qGduyZ6YGp6a6p9/Px6Md0VX275l09PfXpur5l7o6IiCRXKt8BREQkv1QIREQSToVARCThVAhERBJOhUBEJOFUCEREEk6FQCRHFvixme02s8fzlMHN7HXh8++Z2afylONaM3s4H79bep8KQT9hZu8xs2Vm1mRmL5nZ78zsvHzn6mfOAy4C6t39rHyHcffr3f1zxzsfM5ttZo29kakQmFmDmV2X7xxxokLQD5jZR4FvAv8G1AJjgf8ELstnrkxmVpzvDL1gHLDR3fcf74zCrQv9/0lhcHc9YvwAaoAm4Iou2pQRFIot4eObQFk4bTbQCHwM2Aa8BPxtOO1s4GWgKGNe7wRWhM9TwAJgPbATuAMYEk4bDzjwfuAF4KFw/N8Am8L2nwI2Ahf2YH7vC+e3A/hkRq4i4F/C1+4DlgNjwmmTgT8Au4C1wF918V7VAYvCtuuAD4Tj3w8cAtrC9/uzWV5bBHw9zPY8MD/MXBxObwC+APwJOAi8DvhbYE2YeQPw9x3m+Ynwb7IF+Ltwfq8Lp/0E+HxG27cDTwF7gD8Dp2VM2wh8HFgB7AV+CZQDVWGW9nC5moC6LMs2NHxfXgEeBz4HPJwxvdP3GHgrsDpcxheBj2dMuyzM/Er4t5ub8bn+YbjsLwKfJ/wcAtcCDwNfA3aH7/XF4bQvhH+jQ+GyfCff/6NxeOQ9gB7H+QeEuUBremXTSZubgUeBEcDwcCXxuXDa7PD1NwMl4T/tAWBwOH09cFHGvH4FLAiffzicbz1Bsfk+8Itw2vhwpfWzcGVTAUwJ/znPA0rDf+QWjhaCXOb3g3BepwPNwMnh9E8AK4GTAAunDw1/92aCFW4xcCbBinpqJ+/VgwRbU+XANGA78KZw2rWZK78sr70+XOHVA4OB/+O1heAFYGqYpQR4G3BCmPkN4Xt/ZsbfditwSrgcP6eTQhAu1zZgFkFBeh/Byj9d8DcSrMDrgCEExef6jM9AYzefs4UEhbkqzPNi+r3o7j0mWJmfHz4fnLF8ZxEUpYsIvgSMBiaH0+4O//5VBJ/bxwmLZPh3aAE+EC7rPxAUSst4n6/L9/9mnB55D6DHcf4B4b3Ay920WQ+8NWP4LQS7ONIrgYNkFJJwhXJ2+PzzwI/C5wOA/cC4cHgN4UoyHB4V/oMWc3TFPTFj+k2EK/ZwuBI4zNFCkMv86jOmPw5cFT5fC1yWZdmvBP7YYdz3gU9naTuG4NvkgIxxXwR+Ej6/lq4LwQNkfKMHLuS1heDmbv5WdwMfCp//CPhSxrQT6bwQfJewuGe0Xwu8IXy+EfjrjGlfAb6X8RnotBCEK9sWwpV0OO7fOFoIunyPCYrf3wMDs7T59yy/r5agyFdkjLsaWJLxd1jX4XPkwMiM91mFoAcP7aOMv53AsG72wdcR7I5J2xSOOzIPd2/NGD4AVIfPfw68y8zKgHcBT7h7el7jgLvMbI+Z7SFYkbcR/COnbe6Q48iwux8I86flMr+XO8k5hqDgdTQOmJWeZzjf9wIjs7StA3a5+76McZsIvqnm4lXL1+F51nFmdrGZPWpmu8JsbwWGdTK/zL9hR+OAj3VYzjG8+u/c2XvXneEExbizLN29x+8mWK5NZvagmb0+HN/V36wEeCljft8n2DJ4zbKEnyN6sDzSgQpB/D1CsD/0HV202ULwz5U2NhzXLXdfTfBPfzHwHoLCkLaZYN/soIxHubu/mDmLjOcvEew2AcDMKgh23/Rkfp3ZTLCLJdv4BzvMs9rd/yFL2y3AEDMbkDFuLMFukFy8avkIVnQdHXk/wuL6a4JdZLXuPghYTLCbKD2/zHmM7eJ3bwa+0GE5K939Fznk7q4L4u0Euw87y9Lle+zuS939MoIV+d0Eu5jSr+vsb9YMDMuY30B3n5rDsuSyPNKBCkHMuftegl0ut5jZO8ys0sxKwm+aXwmb/QL4VzMbbmbDwvb/3YNf83PgBuACgmMEad8DvmBm4wDC+Xd1ptKdwCVmdo6ZlQKf5ehK71jml+k24HNmNik8I+c0MxsK3AOcaGbXhO9LiZnNNLOTO87A3TcTHD/5opmVm9lpBAeJb88xwx3Ah8xstJkNAv65m/alBMdCtgOtZnYx8OYO87vWzKaYWSXw6S7m9QPgejObFS5/lZm9rUNR68xWYKiZ1WSb6O5twG+Az4SfrykExyDSOn2PzazUzN5rZjXu3kJwULgtfN0Pgb81szeZWSp83ya7+0vAfcDXzWxgOO0EM3tDDsuSXp6JObYVVAj6BXf/BvBR4F8JViqbCc5YuTts8nlgGcEZIyuBJ8JxufoFwX7kB9x9R8b4bxGcSXKfme0jONA7q4ucq4B/JDjw+BLBWSTbCL799Xh+HXyDYMV5H8HK5ocE+5j3EaxcryL4xv8y8GWCFXA2VxMcj9gC3EWwn/sPOWb4Qfj7VwBPEny7b+Xoiu9Vwmw3hLl3E2xxLcqY/juCM7weIDiD6YHOfrG7LyM4ePqdcF7rCPald8vdnyX4G28Id8XUZWk2n2DXy8sExyZ+3GE5unqPrwE2mtkrBAfU/zp83eMEB5j/neCg8YMc3XL9G4JCuTpcnjsJjhnl4lvA5eGFf/+R42sSLX2UXaTPmVk1wamOk9z9+Xzn6W3hN/zvufu4bhuL5JG2CKRPmdkl4e6FKoJ94ysJzmiJPTOrMLO3mlmxmY0m2JVzV75ziXRHhUD62mUcvbBtEsHpn/1ls9QIjnvsJtg1tIbgeIxIQdOuIRGRhNMWgYhIwsWyI7Bhw4b5+PHju223f/9+qqqqog8UkTjnj3N2iHf+OGeHeOcv5OzLly/f4e7Ds02LZSEYP348y5Yt67ZdQ0MDs2fPjj5QROKcP87ZId7545wd4p2/kLObWadXpmvXkIhIwqkQiIgknAqBiEjCRVoIzGyuma01s3VmtiDL9H83s6fCx3NhL4MiItKHIjtYbGZFwC0EN51oBJaa2aKwN0sA3P0jGe3/ETgjqjwiIpJdlFsEZxHcPGKDux8m6Gisq54krybo+EpERPpQZFcWm9nlBPcfvS4cvgaY5e7zs7QdR3iLwrDL22zzmwfMA6itrZ2+cOHCbjM0NTVRXR3fe1XEOX+cs0O888c5O8Q7fyFnnzNnznJ3n5FtWpTXEViWcZ1VnauAOzsrAgDufitwK8CMGTM8l3N1M8/pbWt3vtuwjtPqB3HBiVmvqSg4hXxOcnfinB3inT/O2SHe+eOaPcpdQ428+o5G9XR+V6yriHi3UFHKuPWhDfzfmq1R/hoRkdiJshAsBSaZ2YTwblRXkXHTjTQzOwkYTHDLxUiNHlzJi7sPRv1rRERiJbJCEN4MfT5wL0F3vHe4+yozu9nMLs1oejWwsC+6Ih49qIIX96gQiIhkirSvIXdfTHC7vsxxN3UY/kyUGTLVD67g0Q07cXfMsh3CEBFJnkRdWVw/uIKm5lZeOdia7ygiIgUjUYVg9KAKABr3HMhzEhGRwpGsQjA4LAQ6YCwickSyCkG4RaAzh0REjkpUIRhSVUpFSZHOHBIRyZCoQmBmjB5coS0CEZEMiSoEEOwe0sFiEZGjklcItEUgIvIqiSsE9YMr2H2ghf3NupZARAQSWAiOnDmkA8YiIkACC0H9YJ1CKiKSKXGFYPSgSgAatUUgIgIksBCMGFBGSZFpi0BEJJS4QpBKGXWDKmjcrVNIRUQggYUAdF8CEZFMyS0E2jUkIgIktBDUD65k275mmlvb8h1FRCTvElkI0t1Rb9lzKM9JRETyL5mFQN1Ri4gckchCcOSiMnU+JyISbSEws7lmttbM1pnZgk7a/JWZrTazVWb28yjzpI2sKSdlulOZiAhAcVQzNrMi4BbgIqARWGpmi9x9dUabScCNwLnuvtvMRkSVJ1NJUYqRA8u1a0hEhGi3CM4C1rn7Bnc/DCwELuvQ5gPALe6+G8Ddt0WY51VGD65QNxMiIoC5ezQzNrscmOvu14XD1wCz3H1+Rpu7geeAc4Ei4DPu/vtO5jcPmAdQW1s7feHChd1maGpqorq6Ouu07684xHO72vn67MoeLVdf6ip/oYtzdoh3/jhnh3jnL+Tsc+bMWe7uM7JNi2zXEGBZxnWsOsXAJGA2UA/80cxOcfc9r3mh+63ArQAzZszw2bNndxugoaGBztota17L4w+u57zzL6C4qDCPmXeVv9DFOTvEO3+cs0O888c1e5RrwEZgTMZwPbAlS5vfunuLuz8PrCUoDJEbPbiCtnbn5Vd0LYGIJFuUhWApMMnMJphZKXAVsKhDm7uBOQBmNgw4EdgQYaYjdC2BiEggskLg7q3AfOBeYA1wh7uvMrObzezSsNm9wE4zWw0sAT7h7jujypTp6LUEKgQikmxRHiPA3RcDizuMuynjuQMfDR99qi7cItC1BCKSdIV5lLQPlJcUMay6TLuGRCTxElsIIDhgrF1DIpJ0iS4E9YN1pzIRkWQXgkEVbNlziPb2aC6qExGJg0QXgtGDKzjc1s6OpuZ8RxERyZtkF4L0mUM6TiAiCZboQjB2SNDP0As7dZxARJIr2YVgaCUpgw3bm/IdRUQkbxJdCMqKi6gfXMn6HfvzHUVEJG8SXQgAJg6vYsN2FQIRSS4VgmHVPL+jSaeQikhiqRAMr+JQSzsvqTtqEUkoFYLhVYAOGItIciW+EJwwPLitnI4TiEhSJb4QjBhQRlVpkbYIRCSxEl8IzIyJw6vZoFNIRSShEl8IQKeQikiyqRAQnEL64p6DHDzclu8oIiJ9ToWAo2cOPa/dQyKSQJEWAjOba2ZrzWydmS3IMv1aM9tuZk+Fj+uizNOZI6eQ7tABYxFJnshuXm9mRcAtwEVAI7DUzBa5++oOTX/p7vOjypGLCcPS1xJoi0BEkifKLYKzgHXuvsHdDwMLgcsi/H3HrLK0mLqacp1CKiKJZO7R9LFjZpcDc939unD4GmBW5rd/M7sW+CKwHXgO+Ii7b+5kfvOAeQC1tbXTFy5c2G2GpqYmqqurc8r71aUHOdACnz6nIqf2faEn+QtNnLNDvPPHOTvEO38hZ58zZ85yd5+RdaK7R/IArgBuyxi+Bvh2hzZDgbLw+fXAA7nMe/r06Z6LJUuW5NTO3f1Td6/0qTf93tvb23N+TdR6kr/QxDm7e7zzxzm7e7zzF3J2YJl3sk6NctdQIzAmY7ge2NKhCO109/QNg38ATI8wT5cmDquiqbmV7ft0/2IRSZYoC8FSYJKZTTCzUuAqYFFmAzMblTF4KbAmwjxdmhj2ObReB4xFJGEiKwTu3grMB+4lWMHf4e6rzOxmM7s0bHaDma0ys6eBG4Bro8rTHZ1CKiJJFdnpowDuvhhY3GHcTRnPbwRujDJDrupqKigvSekUUhFJHF1ZHEqljPFDq3QKqYgkjgpBhhPUC6mIJJAKQYYThlexedcBmlvV+ZyIJIcKQYaJw6tpd3hh54F8RxER6TMqBBnSZw7pFFIRSRIVggxHOp/TKaQikiAqBBkGlJcwYkCZTiEVkURRIegguG2ltghEJDlUCDrQjexFJGlUCDo4YXg1ew60sKNJnc+JSDKoEHRw8qgBAKze8kqek4iI9A0Vgg6mjqoBYJUKgYgkhApBBzWVJYweVMHql1QIRCQZVAiymFo3kFVb9uY7hohIn1AhyGJqXQ3P79jP/ubWfEcREYlct/cjMLNPA7nc4b7B3R86/kj5N6VuIO7w7Mv7mD5ucL7jiIhEKpcb02zMcV57jiNHQZlaNxCA1Vv2qhCISL/XbSFw95/2RZBCMqqmnMGVJTpzSEQSIZddQ2NznNced+8Xa04zY0rdQJ05JCKJkMuuoZ8SHCOwLto48BPgZ72QqSBMravhJ3/eSEtbOyVFOqYuIv1XLruG5hzrzM1sLvAtoAi4zd2/1Em7y4FfATPdfdmx/r7eNGXUQA63trNh+35OGjkg33FERCKT81ddM/uvXMZlTCsCbgEuBqYAV5vZlCztBgA3AI/lmqUvpA8Y63oCEenverLPY2rmgJkVA9O7aH8WsM7dN7j7YWAhcFmWdp8DvgIc6kGWyE0cXk15SUoHjEWk38vlYPGNwL8AFWaWXisacBi4tYuXjgY2Zww3ArM6zPsMYIy732NmH+8mxzxgHkBtbS0NDQ3dRaepqSmndp2pq4Q/r95EQ/W2Y57H8Tje/PkU5+wQ7/xxzg7xzh/b7O6e0wP4Yq5tw/ZXEBwXSA9fA3w7YzgFNADjw+EGYEYu854+fbrnYsmSJTm168yNv1nhp376997e3n5c8zlWx5s/n+Kc3T3e+eOc3T3e+Qs5O7DMO1mn5rxryN1vNLPRZnaOmV2QfnTxkkZgTMZwPbAlY3gAcArQYGYbgbOBRWY2I9dMUZtaN5BXDrXSuPtgvqOIiEQml9NHATCzLwFXAauBtnC0A511K7EUmGRmE4AXw9e+Jz3R3fcCwzLm3wB83AvkrCEIzhwCWP3SK4wZUpnnNCIi0ci5EADvBE5y95xu3eXurWY2H7iX4PTRH7n7KjO7mWATZVHP4/atySMHkrLg3gRvmToy33FERCLRk0KwASgBcr6Ho7svBhZ3GHdTJ21n9yBLn6goLeKE4dW6W5mI9Gs9KQQHgKfM7H4yioG739DrqQrIlLqBLH1+V75jiIhEpieFYFH4SJSpdQP57VNb2L3/MIOrSvMdR0Sk1+VcCNz9p2ZWAYx197URZiooU+uCexivfukVzn3dsG5ai4jET0+6mLgEeAr4fTg8zcz6/RZC+swhdTUhIv1VT7qY+AxBtxF7ANz9KWBCBJkKyuCqUupqytXVhIj0Wz0pBK3huf+ZcrmFZexNqathZaO2CESkf+pJIXjGzN4DFJnZJDP7NvDniHIVlBnjB7Nhx352NOV85qyISGz0pBD8I0EPpM3AL4BXgA9HEarQzBwf3Ld42cbdeU4iItL7etLX0AF3/6S7z3T3GeHzguo6Oiqnjh5EWXGKpRt1PYGI9D+5dEP9TXf/sJn9D1mOCbj7pZEkKyClxSmmjRnEMhUCEemHcrmOIH0Xsq9FGaTQzRw/hO8+uJ79za1UlfXkOjwRkcKWyz2Ll4c/H4w+TuGaOWEI31myjidf2MN5k3RhmYj0H7nsGlpJF6eJuvtpvZqoQJ05dhApg6Ubd6kQiEi/kss+jreHPz8Y/kzvKnovQUd0iTCgvISTRw3UAWMR6Xe6PWvI3Te5+ybgXHf/J3dfGT4WAG+JPmLhmDl+CE++sIeWtvZ8RxER6TU9uY6gyszOSw+Y2TlAVe9HKlwzxw/hYEubupsQkX6lJ6e/vB/4kZnVhMN7gL/r/UiFK31h2dLndzFtzKA8pxER6R09uaBsubufDpwGnO7u09z9ieiiFZ4RA8sZN7SSx3WcQET6kR6dEG9mbyPoZqLczABw95sjyFWwZo4fwv1rtuLupN8DEZE468n9CL4HXEnQ55ABVwDjunnNXDNba2brzGxBlunXm9lKM3vKzB42syk9zN/nZo4fzO4DLazf3pTvKCIivaInB4vPcfe/AXa7+2eB1wNjOmtsZkXALcDFwBTg6iwr+p+7+6nuPg34CvCNHqXPg5njhwDw+PPqgE5E+oeeFIJ0B3MHzKwOaKHrG9OcBaxz9w3ufhhYCFyW2cDdM0+/qSIG9zeYMKyKYdWl6ndIRPqNnhwj+B8zGwR8FXiCYKX9gy7ajwY2Zww3ArM6NjKzDwIfBUqBN3Y2MzObB8wDqK2tpaGhodvATU1NObXrqXFVbTz07BYaGvb0+rwzRZW/L8Q5O8Q7f5yzQ7zzxza7u3f7INhyOCdjuAyo6eY1VwC3ZQxfA3y7i/bvAX6aS57p06d7LpYsWZJTu5667Y8bfNw/3+Nb9hyIZP5pUeXvC3HO7h7v/HHO7h7v/IWcHVjmnaxTc9o15O7twNczhpv9tbet7KiRVx9DqAe2dNF+IfCOXPLk25HrCXSjGhHpB3pyjOA+M3u35X7O5FJgkplNMLNS4CpgUWYDM5uUMfg24C89yJM3U0YNpKq0iEfW78x3FBGR49aTYwQfJTig22pmhwhOIXV3H5itsbu3mtl84F6gCPiRu68ys5sJNlEWAfPN7EKCA8+7gfcdx7L0meKiFOe+bhgNa7fpegIRib2cC4G7D+jpzN19MbC4w7ibMp5/qKfzLBRvOnkE963eyrMv7+PkUVlroYhILORcCMzszCyj9wKb3L219yLFw5yTRgDwwLPbVAhEJNZ6cozgP4FHCU4Z/UH4fCHwnJm9OYJsBW3EwHJOHV3D/Wu25juKiMhx6Ukh2Aic4e7T3X06MA14BriQ4KrgxHnj5BE8uXkPu/YfzncUEZFj1pNCMNndV6UH3H01QWHY0Pux4uGNk0fgDg1rt+U7iojIMetJIVhrZt81szeEj/8k2C1URnDWT+KcOrqGYdVlPPCsCoGIxFdPCsG1wDrgw8BHgA3huBZgTm8Hi4NUyphz0nAefG67bl8pIrHVkxvTHHT3r7v7O939He7+NXc/4O7t7p7YPpnfdPII9h1qZfkmXWUsIvHU7emjZvZpcusVtMHdHzr+SPFy3qThlBQZDzy7jbMnDs13HBGRHsvlOoKNOc4r2q44C1R1WTGzJgzlgWe38S9vPTnfcUREeiyXQrAkx3klshAAzJk8gs/ds5oXdh5g7NDKfMcREemRXArBTwl2DXXVoY4DPwF+1guZYudNYSF44NmtXHtuV/fqEREpPN0WAndP5BlBPTF+WBUTh1dx/7PbVAhEJHZ6cvqodOGNJ43gsQ272N+cuG6XRCTmVAh6yRtPHsHhtnYefG57vqOIiPSICkEvOWv8EIYPKOPuJ1/MdxQRkR5RIeglxUUp3jGtjiVrt6kTOhGJFRWCXvSuM+tpaXPuWdHVrZlFRAqLCkEvOnnUQE4eNZBfP6HdQyISHyoEvezdZ47m6c17WLctsd0viUjMqBD0skun1ZEyuOvJxnxHERHJSaSFwMzmmtlaM1tnZguyTP+oma02sxVmdr+ZjYsyT18YMaCcC04czl1PvEh7ey599YmI5FdkhcDMioBbgIuBKcDVZjalQ7MngRnufhpwJ/3klpfvOrOeLXsP8ejzO/MdRUSkW1FuEZwFrHP3De5+mOBG95dlNnD3Je5+IBx8FKiPME+fefOUWgaUFfMbHTQWkRgw92h2X5jZ5cBcd78uHL4GmOXu8ztp/x3gZXf/fCfT5wHzAGpra6cvXLiw2wxNTU1UV1cf4xIcnx8908zjL7XyrTmVlBV31V9f5/KZ/3jFOTvEO3+cs0O88xdy9jlz5ix39xnZpuXS++ixyrb2y1p1zOyvgRnAGzqbmbvfCtwKMGPGDJ89e3a3ARoaGsilXRQqxu7koVsf5eDQE3nLGaOPaR75zH+84pwd4p0/ztkh3vnjmj3KXUONwJiM4XrgNVdamdmFwCeBS929OcI8fWrm+CGMGVLBr5/Q2UMiUtiiLARLgUlmNsHMSoGrgEWZDczsDOD7BEVgW4RZ+lwqZbzzjHoeXreDzbsOdP8CEZE8iawQuHsrMB+4F1gD3OHuq8zsZjO7NGz2VaAa+JWZPWVmizqZXSy956yxFKeMHz78fL6jiIh0KspjBLj7YmBxh3E3ZTy/MMrfn28ja8q59PTR/HLpZj70pkkMrirNdyQRkdfQlcURm3fBRA62tPHfj27KdxQRkaxUCCJ20sgBzD5pOD99ZCOHWtryHUdE5DVUCPrAvAsmsqPpsC4wE5GCpELQB14/cSinjq7htj9uUP9DIlJwVAj6gJkx74KJbNixnz+s2ZrvOCIir6JC0EcuPmUkY4ZUcOtDG/IdRUTkVVQI+khxUYrrzpvI8k27Wb5pV77jiIgcoULQh66YUc+gyhK+26CtAhEpHCoEfaiytJj3nzuB/1uzlcc26F4FIlIYVAj62HXnT2RUTTmf+9/VOoNIRAqCCkEfqygtYsHFk3nmxVfUM6mIFAQVgjy49PQ6zhg7iK/cu5b9za35jiMiCadCkAdmxqfePoXt+5r5bsP6fMcRkYRTIciTM8cO5rJpdfzgjxto3K37FYhI/qgQ5NE/z52MGXz592vzHUVEEkyFII/qBlUw74IT+J+nt7Bsoy4yE5H8UCHIs+vfMJGRA8tZ8JuVHDysbqpFpO+pEORZZWkxX73iNNZta+KLv1uT7zgikkAqBAXg/EnDef95E/jZI5tY8uy2fMcRkYSJtBCY2VwzW2tm68xsQZbpF5jZE2bWamaXR5ml0H3iLScxeeQAPnHn0+xoas53HBFJkMgKgZkVAbcAFwNTgKvNbEqHZi8A1wI/jypHXJSXFPGtq87glUOt/NOdK3BX9xMi0jei3CI4C1jn7hvc/TCwELgss4G7b3T3FUB7hDli46SRA1gwdzIPPLuN/37shXzHEZGEsKi+eYa7eua6+3Xh8DXALHefn6XtT4B73P3OLuY3D5gHUFtbO33hwoXdZmhqaqK6uvrYFiBP2t35xrJmntvdxsdPd06sjVf+tDi+95ninD/O2SHe+Qs5+5w5c5a7+4xs04oj/L2WZdwxVx13vxW4FWDGjBk+e/bsbl/T0NBALu0KzSnTD/G2bz/M99Yc5ncXn03twPJ8R+qxuL73aXHOH+fsEO/8cc0e5a6hRmBMxnA9sCXC39dvjBhYzo+vncmBFufaHy9l36GWfEcSkX4sykKwFJhkZhPMrBS4ClgU4e/rV04ZXcMHp5Xx3NZ9/L/bn6ClTYdRRCQakRUCd28F5gP3AmuAO9x9lZndbGaXApjZTDNrBK4Avm9mq6LKE0enDi/mi+88lT/+ZQc3/malziQSkUhEeYwAd18MLO4w7qaM50sJdhlJJ/5q5ttMRMMAAA1PSURBVBhe3HOQb93/F+oGVfDRi07MdyQR6WciLQTSOz584SS27DnIf9z/FwA+cuEkzLIdixcR6TkVghgwM774rlMxg/+4/y/s2t/MZy89haKUioGIHD8VgpgoLkrx5XefxpCqMr734Hp272/hG1eeTllxUb6jiUjMqRDEiJmx4OLJDK0q5QuL17Dn4GG+f80Mqsv0ZxSRY6feR2PoAxdM5OtXnM6jG3ZxxfceYcP2pnxHEpEYUyGIqXdPr+eH75vBS3sPcsm3H+auJxvzHUlEYkqFIMZmnzSC333ofKbW1fCRXz7Nx+54mv3NrfmOJSIxo0IQc6NqKvj5B2Zxw5sm8ZsnG7nkOw+zsnFvvmOJSIyoEPQDxUUpPnrRidx+3Sz2N7dy6S0P88m7VrJ7/+F8RxORGFAh6EfOOWEY933kDVx7zngWLt3MnK83cPtjm2hrV9cUItI5FYJ+pqaihE9fMpXFN5zPSbUD+ORdz3DZLQ/TsHab+ioSkaxUCPqpk0YOYOG8s/n21Wewq+kw1/54KZd852F+t/Il2rWFICIZdCVSP2ZmXHJ6HW+ZOpK7nmzkuw3r+Yfbn+CE4VVc/4YTuOT0OspLdGWySNJpiyABSotTXDlzLPd/bDbfvvoMSopSfOLOFZz1hf/jU3c/w8rGvdptJJJg2iJIkKJUsIXw9tNG8cj6ndyxbDN3LNvMfz26ickjB3D59HrePGUkY4dW5juqiPQhFYIEMjPOed0wznndMD57sIVFT2/hjqWb+fz/ruHz/7uGySMH8OYptbx56kim1g1Ul9ci/ZwKQcLVVJRwzdnjuObscWzauZ8/rN7Kfau38p0l6/iPB9YxrLqUWROH8vqJQzl74lBOGF6lwiDSz6gQyBHjhlZx3fkTue78iexsauaBZ7fx5/U7eWT9Tv53xUsADKsuY9qYQZxWX8Op9TWcOrqGYdVleU4uIsdDhUCyGlpdxhUzxnDFjDG4Oy/sOsCjG3by2IZdPN24h/uf3Ur6+PKomnIm1Q5g0ojq4FFbTdNhx9219SASAyoE0i0zY9zQKsYNreLKmWMBaGpuZdWLe1n54l6eeXEv67Y3cftjOznU0n7kdTf+6T7qh1QydkgFY4dUUjeogpEDy6mtKWdUTTnDq8soLtKJayL5FmkhMLO5wLeAIuA2d/9Sh+llwM+A6cBO4Ep33xhlJukd1WXFzJo4lFkThx4Z197uvLjnIH/Zto/7HllB2ZBRvLDrAOu376dh7XaaW9tfNQ8zGFJZyrDqMoZWlzK0uoyhVaXUVJQwqDJ8VJQysKKYAeUlVJcVM6C8mKrSYlK6TadIr4msEJhZEXALcBHQCCw1s0Xuvjqj2fuB3e7+OjO7CvgycGVUmSRaqZQxZkglY4ZUknq5hNmzTzkyrb3d2XXgMC/vPcTWVw7x8iuH2Lr3ENubDrOzqZmd+w+zsnEPO5sOs6+brrTNoLKkiMqyYqpKi6goTf8soqw4+FlenKK8pIiy4hRlJSlKi4rCnylKilOUFhklRSlKi1MUp1KUFBnFRSlKUsHPdXvaGNK4h6KUUZxKUZSColSKIjNSKShOpUilIGUWjjNSFpyimzLD7Og0M7SLTApalFsEZwHr3H0DgJktBC4DMgvBZcBnwud3At8xM3Nd3dTvpFLGsOoyhlWXccromi7btra1s+9QK3sOtrDnwGH2HmyhqbmVpkOt7DvUyr5DLew/3MaBw60cONzG/ubgeVNzKzuaDtPc0sahljYOtrRxuLWd5tZ2Wo+lW41H/3SMS5tdKiwO6cJgZAxzdByWfbylJ4bP06UlXWOM4DXNzc2UP3J/OO21BehI+4xJljHf17TPsiwd55u1zOVY+zo2O3DgAJXLG7p/XQEW1wP7D1D5xIPdtjvW5IMrS7nj+tcf46s7F2UhGA1szhhuBGZ11sbdW81sLzAU2NFxZmY2D5gHUFtbS0NDQ7cBmpqacmpXqOKcv7ezV4ePkQCl4aNLReEj0O5OSzu0tEGrO23t0Jp+uNPm0NbOkZ/7Dx6ktKycNod258jPdvfw59GHO7QTPseDYQfn6HQnfGQbJhiRHpfOS3rY/VXT0m3JHBfOC6ClrJ2Sktas0zJmkfHcXzsyS7vOZPvadjzf5Goq2ikuOtT17zyO+UdpQEU7xamDXbY5nq+55a2HIlknRFkIshW9jm9BLm2Cke63ArcCzJgxw2fPnt1tgIaGBnJpV6jinD/O2SHe+eOcHeKdP67ZozxloxEYkzFcD2zprI2ZFQM1wK4IM4mISAdRFoKlwCQzm2BmpcBVwKIObRYB7wufXw48oOMDIiJ9K7JdQ+E+//nAvQQ7a3/k7qvM7GZgmbsvAn4I/JeZrSPYErgqqjwiIpJdpNcRuPtiYHGHcTdlPD8EXBFlBhER6Zou6xQRSTgVAhGRhFMhEBFJOBUCEZGEszierWlm24FNOTQdRparlGMkzvnjnB3inT/O2SHe+Qs5+zh3H55tQiwLQa7MbJm7z8h3jmMV5/xxzg7xzh/n7BDv/HHNrl1DIiIJp0IgIpJw/b0Q3JrvAMcpzvnjnB3inT/O2SHe+WOZvV8fIxARke719y0CERHphgqBiEjC9dtCYGZzzWytma0zswX5ztMVM/uRmW0zs2cyxg0xsz+Y2V/Cn4PzmbErZjbGzJaY2RozW2VmHwrHF/wymFm5mT1uZk+H2T8bjp9gZo+F2X8ZdqVekMysyMyeNLN7wuE4Zd9oZivN7CkzWxaOK/jPTZqZDTKzO83s2fDz//o45U/rl4XAzIqAW4CLgSnA1WY2Jb+puvQTYG6HcQuA+919EnB/OFyoWoGPufvJwNnAB8P3Ow7L0Ay80d1PB6YBc83sbODLwL+H2XcD789jxu58CFiTMRyn7ABz3H1axvn3cfjcpH0L+L27TwZOJ/g7xCl/wN373QN4PXBvxvCNwI35ztVN5vHAMxnDa4FR4fNRwNp8Z+zBsvwWuChuywBUAk8Q3Ft7B1Cc7fNUSA+CO//dD7wRuIfg9q+xyB7m2wgM6zAuFp8bYCDwPOFJN3HLn/nol1sEwGhgc8ZwYzguTmrd/SWA8OeIPOfJiZmNB84AHiMmyxDuWnkK2Ab8AVgP7HH31rBJIX9+vgn8E9AeDg8lPtkhuEf5fWa23MzmheNi8bkBJgLbgR+Hu+ZuM7Mq4pP/iP5aCCzLOJ0nGzEzqwZ+DXzY3V/Jd55cuXubu08j+HZ9FnBytmZ9m6p7ZvZ2YJu7L88cnaVpwWXPcK67n0mwG/eDZnZBvgP1QDFwJvBddz8D2E8cdgNl0V8LQSMwJmO4HtiSpyzHaquZjQIIf27Lc54umVkJQRG43d1/E46O1TK4+x6ggeA4xyAzS9/Br1A/P+cCl5rZRmAhwe6hbxKP7AC4+5bw5zbgLoJCHJfPTSPQ6O6PhcN3EhSGuOQ/or8WgqXApPDsiVKCeyEvynOmnloEvC98/j6C/e4FycyM4P7Ta9z9GxmTCn4ZzGy4mQ0Kn1cAFxIc8FsCXB42K8js7n6ju9e7+3iCz/gD7v5eYpAdwMyqzGxA+jnwZuAZYvC5AXD3l4HNZnZSOOpNwGpikv9V8n2QIsIDOW8FniPY3/vJfOfpJusvgJeAFoJvGe8n2Nd7P/CX8OeQfOfsIv95BLsfVgBPhY+3xmEZgNOAJ8PszwA3heMnAo8D64BfAWX5ztrNcswG7olT9jDn0+FjVfr/NA6fm4xlmAYsCz8/dwOD45Q//VAXEyIiCddfdw2JiEiOVAhERBJOhUBEJOFUCEREEk6FQEQk4VQIJHHMrCn8Od7M3tPL8/6XDsN/7s35i0RBhUCSbDzQo0IQ9mzblVcVAnc/p4eZRPqcCoEk2ZeA88O+8D8Sdj73VTNbamYrzOzvAcxsdni/hZ8DK8Nxd4cdpa1Kd5ZmZl8CKsL53R6OS299WDjvZ8L+96/MmHdDRp/2t4dXamNmXzKz1WGWr/X5uyOJUdx9E5F+awHwcXd/O0C4Qt/r7jPNrAz4k5ndF7Y9CzjF3Z8Ph//O3XeF3VIsNbNfu/sCM5vvQQd2Hb2L4CrU04Fh4WseCqedAUwl6BPoT8C5ZrYaeCcw2d093Q2GSBS0RSBy1JuBvwm7pH6MoKuASeG0xzOKAMANZvY08ChBB4eT6Np5wC886Ol0K/AgMDNj3o3u3k7QPcd44BXgEHCbmb0LOHDcSyfSCRUCkaMM+EcP7pY1zd0nuHt6i2D/kUZmswk6p3u9B3c2exIoz2HenWnOeN5GcFOZVoKtkF8D7wB+36MlEekBFQJJsn3AgIzhe4F/CLvUxsxODHvF7KgG2O3uB8xsMkG31Wkt6dd38BBwZXgcYjhwAUHHcFmF93aocffFwIcJdiuJRELHCCTJVgCt4S6enxDcf3Y88ER4wHY7wbfxjn4PXG9mKwhuS/hoxrRbgRVm9oQHXUKn3UVw28inCXpq/Sd3fzksJNkMAH5rZuUEWxMfObZFFOmeeh8VEUk47RoSEUk4FQIRkYRTIRARSTgVAhGRhFMhEBFJOBUCEZGEUyEQEUm4/w/maKFD0AqLdgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "theta_est = VanillaLinearRegressionGeneralized(X,yNoisy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.61444751, 0.76669669, 0.01412282])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta_est"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Performance: RMSE_start: 1.9583195479394473, RMSE_model: 2.177504859184213\n"
     ]
    }
   ],
   "source": [
    "y_pred = X.dot(theta_est[:m])+theta_est[m]\n",
    "rmse_original = np.sqrt(((yTrue-yNoisy).T).dot(yTrue-yNoisy)/Npoints)\n",
    "rmse_pred = np.sqrt(((yNoisy-y_pred).T).dot(yNoisy-y_pred)/Npoints)\n",
    "print(f'Performance: RMSE_start: {rmse_original}, RMSE_model: {rmse_pred}')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.61444755 0.76669653 0.01412242]\n"
     ]
    }
   ],
   "source": [
    "linReg = LinearRegression()\n",
    "linReg.fit(X,yNoisy)\n",
    "theta_skl = np.concatenate((np.array([linReg.intercept_]),linReg.coef_))\n",
    "print(theta_skl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}