{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f33be53b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x205cf9765b0>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 16184 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "x = np.linspace(-5,6,100)\n",
    "np.random.shuffle(x)\n",
    "\n",
    "y = 5 + x - 2*x**2 + 1e-1*x**3 + 1e-2*x**5\n",
    "y = y + 1e1*np.random.randn(100)\n",
    "\n",
    "x = x.reshape(-1,1)\n",
    "y = y.reshape(-1,1)\n",
    "\n",
    "plt.figure()\n",
    "plt.scatter(x,y,s = 10, c = 'gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2ddbdc10",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x205d23f0ee0>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 18964 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 19896 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20200 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20684 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20892 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20720 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20684 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20616 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20728 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20956 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 21216 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 21552 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20680 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20688 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 21096 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20928 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20824 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20944 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20972 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 37568 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N,D = x.shape\n",
    "\n",
    "sse = np.zeros(20)\n",
    "sse[:] = np.nan\n",
    "\n",
    "theta_norm2 = np.zeros(20)\n",
    "theta_norm2[:] = np.nan\n",
    "\n",
    "# feature 변환을 한 이후의 design matrix\n",
    "Z = np.ones((N,1))\n",
    "for D in range(1, 20):\n",
    "    Z = np.hstack([Z, x**D])\n",
    "    \n",
    "    theta = np.linalg.inv(Z.T.dot(Z)).dot(Z.T).dot(y)\n",
    "    y_pred = Z.dot(theta)\n",
    "    \n",
    "    sse[D] = np.sum( (y - y_pred)**2 )    \n",
    "    theta_norm2[D] = np.sum( theta**2 )\n",
    "    \n",
    "    plt.figure()\n",
    "    plt.scatter(x,y,s = 10, c = 'gray')\n",
    "    plt.scatter(x,y_pred,s = 10, c = 'red')\n",
    "    \n",
    "plt.figure()\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(sse)\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(theta_norm2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "3e6fa20e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x205d20a9d90>]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 19012 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 19960 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20500 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20448 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20732 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20348 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20316 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20708 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20800 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20856 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 21388 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 21112 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20848 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20748 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20928 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20528 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20696 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20716 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20648 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 27904 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N,D = x.shape\n",
    "alpha = 1e2\n",
    "\n",
    "sse = np.zeros(20)\n",
    "sse[:] = np.nan\n",
    "\n",
    "theta_norm2 = np.zeros(20)\n",
    "theta_norm2[:] = np.nan\n",
    "\n",
    "# ridge regression, normal equation\n",
    "# feature 변환을 한 이후의 design matrix\n",
    "Z = np.ones((N,1))\n",
    "for D in range(1, 20):\n",
    "    Z = np.hstack([Z, x**D])\n",
    "    \n",
    "    theta = np.linalg.inv(Z.T.dot(Z) + (alpha/2)*np.identity(D + 1)).dot(Z.T).dot(y)\n",
    "    y_pred = Z.dot(theta)\n",
    "    \n",
    "    sse[D] = np.sum( (y - y_pred)**2 )    \n",
    "    theta_norm2[D] = np.sum( theta**2 )\n",
    "    \n",
    "    plt.figure()\n",
    "    plt.scatter(x,y,s = 10, c = 'gray')\n",
    "    plt.scatter(x,y_pred,s = 10, c = 'red')\n",
    "    \n",
    "plt.figure()\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(sse)\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(theta_norm2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "a7b69c2f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "norm squre: linear 7.57 ridge 4.44 lasso 0.48\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\jacob\\anaconda3\\lib\\site-packages\\sklearn\\linear_model\\_coordinate_descent.py:647: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations, check the scale of the features or consider increasing regularisation. Duality gap: 6.066e+03, tolerance: 8.213e+00\n",
      "  model = cd_fast.enet_coordinate_descent(\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x205d5a2ac10>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 25592 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.linear_model import LinearRegression, Ridge, Lasso\n",
    "\n",
    "ridge_alpha = 1e1\n",
    "lasso_alpha = 1e0\n",
    "\n",
    "linear = LinearRegression()\n",
    "ridge = Ridge(alpha = ridge_alpha)\n",
    "lasso = Lasso(alpha = lasso_alpha, fit_intercept = False)\n",
    "\n",
    "Z = np.ones((N,1))\n",
    "for D in range(1,10):\n",
    "    Z = np.hstack([Z, x**D])\n",
    "\n",
    "linear.fit(Z, y)\n",
    "ridge.fit(Z, y)\n",
    "lasso.fit(Z, y)\n",
    "\n",
    "y_linear = linear.predict(Z)\n",
    "y_rigde = linear.predict(Z)\n",
    "y_lasso = linear.predict(Z)\n",
    "\n",
    "linear_norm2 = np.sum( linear.coef_**2 )\n",
    "ridge_norm2 = np.sum( ridge.coef_**2 )\n",
    "lasso_norm2 = np.sum( lasso.coef_**2 )\n",
    "\n",
    "print('norm squre: linear %.2f ridge %.2f lasso %.2f' % (linear_norm2, ridge_norm2, lasso_norm2))\n",
    "\n",
    "plt.figure()\n",
    "plt.scatter(x,y,s = 10, c = 'gray')\n",
    "plt.scatter(x,y_linear,s = 10)\n",
    "plt.scatter(x,y_rigde,s = 10)\n",
    "plt.scatter(x,y_lasso,s = 10)\n",
    "plt.legend(['data','linear','ridge','lasso'])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
