2주차_Regression_과제1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "HtloJibI_5uB"
   },
   "source": [
    "# 과제 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "poKqKrLb_5uH"
   },
   "source": [
    "## Matrix 구현 \n",
    "1. Linear Regression에서 쓰이는 **정규방정식**을 행렬로 구현하고, sklearn 혹은 OLS 패키지를 통해 구한 **실제 값과 비교**해주세요.\n",
    "2. LSE에서 쓰이는 Loss Function, **MSE**를 행렬로 구현해 출력해주세요."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "MaVCkq3g_5uJ"
   },
   "source": [
    "### Data Load "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "id": "lX6_aNNU_5uK"
   },
   "outputs": [],
   "source": [
    "import pandas as pd \n",
    "data = pd.read_csv(\"data/w2_regression_assignment1.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "id": "PJVRecL-_5uM",
    "outputId": "78c70b21-ad58-4a1f-b5b3-f90a6f55d1c7"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>y</th>\n",
       "      <th>x1</th>\n",
       "      <th>x2</th>\n",
       "      <th>x3</th>\n",
       "      <th>x4</th>\n",
       "      <th>x5</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>10</td>\n",
       "      <td>38.9</td>\n",
       "      <td>64.7</td>\n",
       "      <td>4</td>\n",
       "      <td>868</td>\n",
       "      <td>59.7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>13</td>\n",
       "      <td>41.6</td>\n",
       "      <td>45.3</td>\n",
       "      <td>-4</td>\n",
       "      <td>957</td>\n",
       "      <td>61.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>11</td>\n",
       "      <td>39.7</td>\n",
       "      <td>74.1</td>\n",
       "      <td>8</td>\n",
       "      <td>786</td>\n",
       "      <td>61.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>7</td>\n",
       "      <td>37.3</td>\n",
       "      <td>48.0</td>\n",
       "      <td>19</td>\n",
       "      <td>984</td>\n",
       "      <td>67.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>10</td>\n",
       "      <td>39.5</td>\n",
       "      <td>51.9</td>\n",
       "      <td>6</td>\n",
       "      <td>700</td>\n",
       "      <td>57.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>9</td>\n",
       "      <td>37.4</td>\n",
       "      <td>53.6</td>\n",
       "      <td>-5</td>\n",
       "      <td>1037</td>\n",
       "      <td>58.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>9</td>\n",
       "      <td>35.1</td>\n",
       "      <td>71.4</td>\n",
       "      <td>3</td>\n",
       "      <td>986</td>\n",
       "      <td>58.6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>6</td>\n",
       "      <td>38.8</td>\n",
       "      <td>58.3</td>\n",
       "      <td>6</td>\n",
       "      <td>819</td>\n",
       "      <td>59.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>5</td>\n",
       "      <td>36.6</td>\n",
       "      <td>52.6</td>\n",
       "      <td>-19</td>\n",
       "      <td>791</td>\n",
       "      <td>54.4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    y    x1    x2  x3    x4    x5\n",
       "0  10  38.9  64.7   4   868  59.7\n",
       "1  13  41.6  45.3  -4   957  61.4\n",
       "2  11  39.7  74.1   8   786  61.0\n",
       "3   7  37.3  48.0  19   984  67.5\n",
       "4  10  39.5  51.9   6   700  57.2\n",
       "5   9  37.4  53.6  -5  1037  58.8\n",
       "6   9  35.1  71.4   3   986  58.6\n",
       "7   6  38.8  58.3   6   819  59.2\n",
       "8   5  36.6  52.6 -19   791  54.4"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "id": "Cdf4C5No_5uP",
    "outputId": "28e40fcd-7aed-40be-b449-38b4c9937e5c"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9, 6)"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.shape     # y:(9,1), X:(9,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "id": "rZuPCz_9_5uQ"
   },
   "outputs": [],
   "source": [
    "X = data.drop([\"y\"], axis=1)\n",
    "y = data.y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x1</th>\n",
       "      <th>x2</th>\n",
       "      <th>x3</th>\n",
       "      <th>x4</th>\n",
       "      <th>x5</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>38.9</td>\n",
       "      <td>64.7</td>\n",
       "      <td>4</td>\n",
       "      <td>868</td>\n",
       "      <td>59.7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>41.6</td>\n",
       "      <td>45.3</td>\n",
       "      <td>-4</td>\n",
       "      <td>957</td>\n",
       "      <td>61.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>39.7</td>\n",
       "      <td>74.1</td>\n",
       "      <td>8</td>\n",
       "      <td>786</td>\n",
       "      <td>61.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>37.3</td>\n",
       "      <td>48.0</td>\n",
       "      <td>19</td>\n",
       "      <td>984</td>\n",
       "      <td>67.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>39.5</td>\n",
       "      <td>51.9</td>\n",
       "      <td>6</td>\n",
       "      <td>700</td>\n",
       "      <td>57.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>37.4</td>\n",
       "      <td>53.6</td>\n",
       "      <td>-5</td>\n",
       "      <td>1037</td>\n",
       "      <td>58.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>35.1</td>\n",
       "      <td>71.4</td>\n",
       "      <td>3</td>\n",
       "      <td>986</td>\n",
       "      <td>58.6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>38.8</td>\n",
       "      <td>58.3</td>\n",
       "      <td>6</td>\n",
       "      <td>819</td>\n",
       "      <td>59.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>36.6</td>\n",
       "      <td>52.6</td>\n",
       "      <td>-19</td>\n",
       "      <td>791</td>\n",
       "      <td>54.4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     x1    x2  x3    x4    x5\n",
       "0  38.9  64.7   4   868  59.7\n",
       "1  41.6  45.3  -4   957  61.4\n",
       "2  39.7  74.1   8   786  61.0\n",
       "3  37.3  48.0  19   984  67.5\n",
       "4  39.5  51.9   6   700  57.2\n",
       "5  37.4  53.6  -5  1037  58.8\n",
       "6  35.1  71.4   3   986  58.6\n",
       "7  38.8  58.3   6   819  59.2\n",
       "8  36.6  52.6 -19   791  54.4"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    10\n",
       "1    13\n",
       "2    11\n",
       "3     7\n",
       "4    10\n",
       "5     9\n",
       "6     9\n",
       "7     6\n",
       "8     5\n",
       "Name: y, dtype: int64"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "id": "sI24Ruws_5uQ"
   },
   "outputs": [],
   "source": [
    "# Matrix 계산을 위해 X와 y를 numpy 형태로 바꾸어줍니다. \n",
    "X = X.to_numpy()\n",
    "y = y.to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "id": "oIEx1Gk5_5uR",
    "outputId": "af2ee5f2-5512-404c-91b7-a0afcfdbed01"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  38.9,   64.7,    4. ,  868. ,   59.7],\n",
       "       [  41.6,   45.3,   -4. ,  957. ,   61.4],\n",
       "       [  39.7,   74.1,    8. ,  786. ,   61. ],\n",
       "       [  37.3,   48. ,   19. ,  984. ,   67.5],\n",
       "       [  39.5,   51.9,    6. ,  700. ,   57.2],\n",
       "       [  37.4,   53.6,   -5. , 1037. ,   58.8],\n",
       "       [  35.1,   71.4,    3. ,  986. ,   58.6],\n",
       "       [  38.8,   58.3,    6. ,  819. ,   59.2],\n",
       "       [  36.6,   52.6,  -19. ,  791. ,   54.4]])"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "id": "8xH58zzD_5uS",
    "outputId": "353a1f50-ddc0-4f27-9c69-ec6b1b2f2046"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([10, 13, 11,  7, 10,  9,  9,  6,  5], dtype=int64)"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "F3JfJR95_5uT"
   },
   "source": [
    "### 1. Normal Equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "id": "Fs3H6XKI_5uX"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.linalg import inv "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "id": "oEiAbw86_5uX"
   },
   "outputs": [],
   "source": [
    "def estimate_beta(X, y): \n",
    "    # Beta_hat = (X'X)^(-1)X'y\n",
    "    designX = np.concatenate((np.ones(X.shape[0]).reshape(-1,1), X),axis=1)   # designX: X matrix의 첫번째 열에 '1' 벡터 추가\n",
    "    beta_hat = inv(np.transpose(designX)@designX)@np.transpose(designX)@y     \n",
    "    \n",
    "    return beta_hat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-3.92447368e+01,  1.31232583e+00,  8.53744361e-02,  7.41849897e-02,\n",
       "        1.50018573e-02, -3.42273652e-01])"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# X와 y로 회귀계수 추정\n",
    "beta_hat = estimate_beta(X, y)\n",
    "beta_hat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Kr4KXjl9_5uZ"
   },
   "source": [
    "#### sklearn API로 계산한 값과 비교"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "id": "Ef3Ggqo8_5uZ"
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "id": "6Parv6ui_5uZ"
   },
   "outputs": [],
   "source": [
    "model = LinearRegression()\n",
    "model.fit(X,y)\n",
    "\n",
    "pred = model.predict(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-39.24473678135658\n",
      "[ 1.31232583  0.08537444  0.07418499  0.01500186 -0.34227365]\n"
     ]
    }
   ],
   "source": [
    "print(model.intercept_)   # intercept\n",
    "print(model.coef_)   # 추정된 회귀계수 (intercept 제외)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Csj01bXd_5ua"
   },
   "source": [
    "### 2. MSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "id": "ADYxjyQ9_5ub"
   },
   "outputs": [],
   "source": [
    "def MSE(X, y, beta_hat):\n",
    "    designX = np.concatenate((np.ones(X.shape[0]).reshape(-1,1), X),axis=1)  # designX: X matrix의 첫번째 열에 '1' 벡터 추가\n",
    "    y_pred = designX@beta_hat    # Y_hat = XB_hat; beta_hat을 이용하여 Y prediction\n",
    "    mse = np.sum(np.transpose(y-pred)@(y-pred)) / (X.shape[0] - len(beta_hat))   # MSE = sum((Y - Y_hat)^2) / (n-p)\n",
    "    \n",
    "    return mse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.84679316870574"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MSE(X, y, estimate_beta(X, y))    # X: data, Y: true value, estimate_beta(X,y): predicted beta"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "Regression_과제1",
   "provenance": []
  },
  "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}