addendum/07.04.03.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "m7xmMHOY5YUe"
   },
   "source": [
    "# 「7.4.3 当てはまりの良さの指標の問題点」についての補足\n",
    "\n",
    "p.193--194に次のようなコードがあります．この結果が書籍に掲載されているものと同じにならないという指摘をいただきました．このことについて，補足します．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "kgskdjhp5M7W",
    "outputId": "cbbf2015-552e-4ff7-82fc-3e9efb3d9fdf"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearRegression()"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "7.725798977575843e-07"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "\n",
    "my_data = sm.datasets.get_rdataset('cars', 'datasets').data\n",
    "\n",
    "my_idx = [1, 10, 26, 33, 38, 43]\n",
    "my_sample = my_data.iloc[my_idx, ]\n",
    "X, y = my_sample[['speed']], my_sample['dist']\n",
    "\n",
    "d = 5\n",
    "X5 = PolynomialFeatures(d).fit_transform(X) # Xの「0」乗から5乗の変数\n",
    "\n",
    "my_model = LinearRegression()\n",
    "my_model.fit(X5, y)\n",
    "y_ = my_model.predict(X5)\n",
    "\n",
    "((y - y_)**2).mean()**0.5\n",
    "#> 7.725744805546204e-07 # RMSE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "u4CWXUtv5_w2"
   },
   "source": [
    "## 簡単な説明\n",
    "\n",
    "コメントとして掲載しているのは，本書執筆時のGoogle Colabの結果です．環境によっては，これとは違う結果になることがあります．本書で強調しているように，回帰分析は乱数を使わないので，結果はいつも同じになるはずです．それにもかかわらず，違う結果が出るのはおかしい，というわけです．\n",
    "\n",
    "ごもっともです．\n",
    "\n",
    "違う結果が出る原因をひとことで言えば，「浮動小数点数の計算でよくある誤差」です．正しく計算すれば同じ結果になるはずのものでも，コンピュータで計算すると，違う結果になることがあるのです．\n",
    "\n",
    "たとえば，上で作ったX5を左右反転しても結果は変わらないはずですが，数値計算の結果は（おそらく）変わります．（正確な計算の結果は後で示します．）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "qAXpIdEBvmb6",
    "outputId": "c04b0b08-fdee-490f-954d-8ce4e22e8a3d"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearRegression()"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "8.659166545298886e-11"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X5a = np.fliplr(X5) # 左右反転\n",
    "my_model = LinearRegression()\n",
    "my_model.fit(X5a, y)\n",
    "y_ = my_model.predict(X5a)\n",
    "((y - y_)**2).mean()**0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "0cV_jaAIyGb6"
   },
   "source": [
    "以上で納得できる方にとっては，この話はこれで終わりでいいでしょう．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "S5K1j9mBlpN5"
   },
   "source": [
    "## 詳しい説明\n",
    "\n",
    "**話をDockerで試しやすいのものに限定します．**\n",
    "\n",
    "scikit-learnでは，回帰分析の主要部の計算に，BLASとLAPACK（以下，まとめてBLAS）という線形代数のライブラリを使います．BLASには複数の実装があります．ここでは次のものを対象にします．インストール方法がそれぞれ2通りあります．\n",
    "\n",
    "- OpenBLAS\n",
    "    - `pip install numpy`\n",
    "    - `apt install libopenblas-dev`\n",
    "- Intel MKL\n",
    "    - `conda install numpy`\n",
    "    - `apt install intel-mkl`\n",
    "\n",
    "次のアーキテクチャを想定します．ただし，aarch64ではIntel MKLは使えません．\n",
    "\n",
    "- x86_64: `docker run -it --rm ubuntu`\n",
    "  - Intel\n",
    "  - AMD\n",
    "- Apple\n",
    "  - Rosetta 2: `docker run -it --rm --platform linux/x86_64 ubuntu`（aarch64によるx86_64のエミュレーション）\n",
    "  - aarch64: `docker run -it --rm --platform linux/arm64 ubuntu` \n",
    "\n",
    "計算結果は，BLASとアーキテクチャの組合せによります．組合せが異なれば，結果も異なると考えてください．次のような，例外的な事例があります．\n",
    "\n",
    "- Intel MKLを使う場合は，[環境変数`MKL_CBWR`](https://jp.xlsoft.com/documents/intel/mkl/11.3/mkl113_userguide_win/GUID-DCB010F6-DDBF-4A00-8BB3-049BEFDC2ED2.htm)に注意が必要でです．\n",
    "  - 値を`COMPATIBLE`に設定すると，Intel, AMD, Rosetta 2で同じ結果になりました．\n",
    "  - Intel（AVX2をサポートするCPU）では，`MKL_CBWR`を設定していない場合は`AUTO`や`AVX2`と同じ（`COMPATIBLE`とは異なる）結果でした．\n",
    "  - AMDとRosetta 2は同じ結果でした．`MKL_CBWR`を設定していない場合は`COMPATIBLE`とは異なる結果，設定しいる場合は`COMPATIBLE`と同じ結果になりました．\n",
    "- x86_64では，OpenBLASのバージョンによって，結果が変わることがありました．aarch64やRosetta 2ではそういうことはなかったので，AVXを使う場合だけのことかもしれません．OpenBLASをnumpyと一緒にインストールする場合は，numpyのバージョンによってOpenBLASのバージョンが変わることがあるので注意が必要です．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 再現方法\n",
    "\n",
    "上述の方法で構築したDockerコンテナでの再現方法をまとめます．\n",
    "補足：numpyとscipyのバージョンは，執筆時点のMinicondaに合わせています．\n",
    "\n",
    "#### OpenBLAS\n",
    "\n",
    "##### numpyと一緒にインストールする場合\n",
    "\n",
    "```\n",
    "apt update\n",
    "apt install -y python3 python3-pip\n",
    "pip install numpy==1.21.2 scipy==1.7.3 scikit-learn statsmodels\n",
    "python3 # Pythonを起動し，計算する．（以下略）\n",
    "```\n",
    "\n",
    "##### numpyとは別にインストールする場合\n",
    "\n",
    "```\n",
    "apt update\n",
    "apt install -y python3 python3-pip libopenblas-dev\n",
    "pip install --no-binary :all: numpy==1.21.2\n",
    "pip install scipy==1.7.3 scikit-learn statsmodels\n",
    "python3 # Pythonを起動し，計算する．（以下略）\n",
    "```\n",
    "\n",
    "#### Intel MKL\n",
    "\n",
    "##### numpyと一緒にインストールする場合\n",
    "\n",
    "```\n",
    "apt update\n",
    "apt install -y curl\n",
    "curl -O https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh\n",
    "sh Miniconda3-latest-Linux-x86_64.sh -b\n",
    "eval \"$(~/miniconda3/bin/conda shell.bash hook)\"\n",
    "conda install -y scikit-learn pandas statsmodels\n",
    "\n",
    "python3 # Pythonを起動し，計算する．（以下略）\n",
    "# あるいは\n",
    "MKL_CBWE=COMPATIBLE python3 # Pythonを起動し，計算する．（以下略）\n",
    "```\n",
    "\n",
    "##### numpyとは別にインストールする場合\n",
    "\n",
    "```\n",
    "apt update\n",
    "apt install -y python3 python3-pip intel-mkl\n",
    "pip install --no-binary :all: numpy==1.21.2\n",
    "pip install scipy==1.7.3 scikit-learn statsmodels\n",
    "\n",
    "python3 # Pythonを起動し，計算する．（以下略）\n",
    "# あるいは\n",
    "MKL_CBWE=COMPATIBLE python3 # Pythonを起動し，計算する．（以下略）\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "pAj1wvpd9F22"
   },
   "source": [
    "## お詫び\n",
    "\n",
    "上のコードの`X5`には少し問題がありました．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 238
    },
    "id": "lqd_10R59ciz",
    "outputId": "c8ff7e56-ab94-4b2c-e7c9-dfcb22c849a7"
   },
   "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>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>64.0</td>\n",
       "      <td>256.0</td>\n",
       "      <td>1024.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0</td>\n",
       "      <td>11.0</td>\n",
       "      <td>121.0</td>\n",
       "      <td>1331.0</td>\n",
       "      <td>14641.0</td>\n",
       "      <td>161051.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>256.0</td>\n",
       "      <td>4096.0</td>\n",
       "      <td>65536.0</td>\n",
       "      <td>1048576.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.0</td>\n",
       "      <td>18.0</td>\n",
       "      <td>324.0</td>\n",
       "      <td>5832.0</td>\n",
       "      <td>104976.0</td>\n",
       "      <td>1889568.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>400.0</td>\n",
       "      <td>8000.0</td>\n",
       "      <td>160000.0</td>\n",
       "      <td>3200000.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1.0</td>\n",
       "      <td>22.0</td>\n",
       "      <td>484.0</td>\n",
       "      <td>10648.0</td>\n",
       "      <td>234256.0</td>\n",
       "      <td>5153632.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     0     1      2        3         4          5\n",
       "0  1.0   4.0   16.0     64.0     256.0     1024.0\n",
       "1  1.0  11.0  121.0   1331.0   14641.0   161051.0\n",
       "2  1.0  16.0  256.0   4096.0   65536.0  1048576.0\n",
       "3  1.0  18.0  324.0   5832.0  104976.0  1889568.0\n",
       "4  1.0  20.0  400.0   8000.0  160000.0  3200000.0\n",
       "5  1.0  22.0  484.0  10648.0  234256.0  5153632.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame(X5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "qJGQ2E3a9gW1"
   },
   "source": [
    "左端のすべてが1の列は，線形モデルの定数に相当する部分です．`LinearRegression()`で改めて切片を導入するので，この列は不要です．\n",
    "\n",
    "この列の有無の予測値への影響は，正確な計算なら無いのですが，数値計算では有ります．（正確な計算の結果は後で示します．）\n",
    "\n",
    "ですから，次のようにすべきでした．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 238
    },
    "id": "z1gdH2qi-sQO",
    "outputId": "b021af44-02b4-4b9c-f7f2-f26fc0cedef3"
   },
   "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>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>4.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>64.0</td>\n",
       "      <td>256.0</td>\n",
       "      <td>1024.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>11.0</td>\n",
       "      <td>121.0</td>\n",
       "      <td>1331.0</td>\n",
       "      <td>14641.0</td>\n",
       "      <td>161051.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>16.0</td>\n",
       "      <td>256.0</td>\n",
       "      <td>4096.0</td>\n",
       "      <td>65536.0</td>\n",
       "      <td>1048576.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>18.0</td>\n",
       "      <td>324.0</td>\n",
       "      <td>5832.0</td>\n",
       "      <td>104976.0</td>\n",
       "      <td>1889568.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>20.0</td>\n",
       "      <td>400.0</td>\n",
       "      <td>8000.0</td>\n",
       "      <td>160000.0</td>\n",
       "      <td>3200000.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>22.0</td>\n",
       "      <td>484.0</td>\n",
       "      <td>10648.0</td>\n",
       "      <td>234256.0</td>\n",
       "      <td>5153632.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      0      1        2         3          4\n",
       "0   4.0   16.0     64.0     256.0     1024.0\n",
       "1  11.0  121.0   1331.0   14641.0   161051.0\n",
       "2  16.0  256.0   4096.0   65536.0  1048576.0\n",
       "3  18.0  324.0   5832.0  104976.0  1889568.0\n",
       "4  20.0  400.0   8000.0  160000.0  3200000.0\n",
       "5  22.0  484.0  10648.0  234256.0  5153632.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X5 = PolynomialFeatures(d, include_bias=False).fit_transform(X) # Xの1乗から5乗の変数\n",
    "pd.DataFrame(X5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fQGjpGMiVA7w"
   },
   "source": [
    "RMSE（訓練）はやはりほぼ0になります．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "WPwGffNP-79E",
    "outputId": "9f16c5c1-0b5e-4f7d-f4ee-6fd6dae094ed"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearRegression()"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "4.943149344979802e-07"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_model = LinearRegression()\n",
    "my_model.fit(X5, y)\n",
    "y_ = my_model.predict(X5)\n",
    "((y - y_)**2).mean()**0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sEe68lGT_j6c"
   },
   "source": [
    "## scikit-learnがやっていること\n",
    "\n",
    "`my_model.fit(X5, y)`と`my_model.predict(X5)`で行われるのは，次のような処理です（先と同じ結果になります）．\n",
    "\n",
    "BLASとアーキテクチャの組合せによって結果が変わるのは★の部分です．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "l9-kv_Eb6IMl",
    "outputId": "9815fda8-1b27-428a-9a02-b53e21a61809"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.943149344979802e-07"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import scipy\n",
    "M = np.insert(X5, 5, 1, axis=1)\n",
    "X_offset = M.mean(axis=0)\n",
    "y_offset = y.mean()\n",
    "b, _, _, _ = scipy.linalg.lstsq(M - X_offset, y - y_offset) # ★\n",
    "y_ = (y_offset -\n",
    "      np.dot(b, X_offset) + # ★\n",
    "      X5 @ b[:5])           # ★\n",
    "((y - y_)**2).mean()**0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "BfxgVxmoYDX9"
   },
   "source": [
    "Mはp.230の脚註4にあるような行列です．整形して表示します．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 238
    },
    "id": "UrpAR7nUYClv",
    "outputId": "3e1ba531-7150-4fdd-84ea-6c78af65b2c4"
   },
   "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>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>4.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>64.0</td>\n",
       "      <td>256.0</td>\n",
       "      <td>1024.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>11.0</td>\n",
       "      <td>121.0</td>\n",
       "      <td>1331.0</td>\n",
       "      <td>14641.0</td>\n",
       "      <td>161051.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>16.0</td>\n",
       "      <td>256.0</td>\n",
       "      <td>4096.0</td>\n",
       "      <td>65536.0</td>\n",
       "      <td>1048576.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>18.0</td>\n",
       "      <td>324.0</td>\n",
       "      <td>5832.0</td>\n",
       "      <td>104976.0</td>\n",
       "      <td>1889568.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>20.0</td>\n",
       "      <td>400.0</td>\n",
       "      <td>8000.0</td>\n",
       "      <td>160000.0</td>\n",
       "      <td>3200000.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>22.0</td>\n",
       "      <td>484.0</td>\n",
       "      <td>10648.0</td>\n",
       "      <td>234256.0</td>\n",
       "      <td>5153632.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      0      1        2         3          4    5\n",
       "0   4.0   16.0     64.0     256.0     1024.0  1.0\n",
       "1  11.0  121.0   1331.0   14641.0   161051.0  1.0\n",
       "2  16.0  256.0   4096.0   65536.0  1048576.0  1.0\n",
       "3  18.0  324.0   5832.0  104976.0  1889568.0  1.0\n",
       "4  20.0  400.0   8000.0  160000.0  3200000.0  1.0\n",
       "5  22.0  484.0  10648.0  234256.0  5153632.0  1.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "82VsR5PWghzU"
   },
   "source": [
    "## 正確な計算\n",
    "\n",
    "正確に計算すると，RMSE（訓練）は0になります．Wolfram|Alphaを使うのが簡単です．\n",
    "\n",
    "[degree 5 {4,10},{11,28},{16,32},{18,76},{20,32},{22,66}](https://www.wolframalpha.com/input/?i=degree+5+%7B4%2C10%7D%2C%7B11%2C28%7D%2C%7B16%2C32%7D%2C%7B18%2C76%7D%2C%7B20%2C32%7D%2C%7B22%2C66%7D)\n",
    "\n",
    "計算方法はいろいろありますが，ここでは，ムーア・ペンローズ型一般逆行列（p.230の脚註4）を使います．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 139
    },
    "id": "OD4mMz6-dbsW",
    "outputId": "06450323-96aa-43f5-b864-0a287070b437"
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}4 & 16 & 64 & 256 & 1024 & 1\\\\11 & 121 & 1331 & 14641 & 161051 & 1\\\\16 & 256 & 4096 & 65536 & 1048576 & 1\\\\18 & 324 & 5832 & 104976 & 1889568 & 1\\\\20 & 400 & 8000 & 160000 & 3200000 & 1\\\\22 & 484 & 10648 & 234256 & 5153632 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[ 4,  16,    64,    256,    1024, 1],\n",
       "[11, 121,  1331,  14641,  161051, 1],\n",
       "[16, 256,  4096,  65536, 1048576, 1],\n",
       "[18, 324,  5832, 104976, 1889568, 1],\n",
       "[20, 400,  8000, 160000, 3200000, 1],\n",
       "[22, 484, 10648, 234256, 5153632, 1]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import *\n",
    "Y = Matrix(y)\n",
    "M = Matrix(np.array(X.speed)[:, np.newaxis]**[1, 2, 3, 4, 5, 0])\n",
    "M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 139
    },
    "id": "xmibIwHffrxx",
    "outputId": "10e25544-0006-4d87-b04b-d0d9ac99eb1d"
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{2786331301}{582120}\\\\- \\frac{126621007}{155232}\\\\\\frac{919150}{14553}\\\\- \\frac{1420541}{620928}\\\\\\frac{293749}{9313920}\\\\- \\frac{1407358}{147}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[2786331301/582120],\n",
       "[-126621007/155232],\n",
       "[     919150/14553],\n",
       "[  -1420541/620928],\n",
       "[   293749/9313920],\n",
       "[     -1407358/147]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = M.pinv() * Y\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 37
    },
    "id": "eKdNIdm2fwHm",
    "outputId": "5069e10c-12df-4280-c7c7-cc090135c7b0"
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y_ = M * b\n",
    "(Y - Y_).norm()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nhL0wQU6h0gH"
   },
   "source": [
    "この方法なら，「お詫び」で書いたように，すべて1の列が2列あっても，予測結果は変わりません．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 139
    },
    "id": "LXPYrwotgNio",
    "outputId": "a0bacad5-fe31-4164-c01d-297e4df2f5b6"
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 4 & 16 & 64 & 256 & 1024 & 1\\\\1 & 11 & 121 & 1331 & 14641 & 161051 & 1\\\\1 & 16 & 256 & 4096 & 65536 & 1048576 & 1\\\\1 & 18 & 324 & 5832 & 104976 & 1889568 & 1\\\\1 & 20 & 400 & 8000 & 160000 & 3200000 & 1\\\\1 & 22 & 484 & 10648 & 234256 & 5153632 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1,  4,  16,    64,    256,    1024, 1],\n",
       "[1, 11, 121,  1331,  14641,  161051, 1],\n",
       "[1, 16, 256,  4096,  65536, 1048576, 1],\n",
       "[1, 18, 324,  5832, 104976, 1889568, 1],\n",
       "[1, 20, 400,  8000, 160000, 3200000, 1],\n",
       "[1, 22, 484, 10648, 234256, 5153632, 1]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = Matrix(np.array(X.speed)[:, np.newaxis]**[0, 1, 2, 3, 4, 5, 0])\n",
    "M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 37
    },
    "id": "yQRSYGxRiNFO",
    "outputId": "0fbb8372-3598-486d-846c-d99444d99f2e"
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = M.pinv() * Y\n",
    "Y_ = M * b\n",
    "(Y - Y_).norm()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "RELPu7etzRFG"
   },
   "source": [
    "左右反転させても同じです．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 139
    },
    "id": "O9-sSPYZyzv_",
    "outputId": "697196cc-548e-4aab-f8f6-61164adfd7b7"
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 1024 & 256 & 64 & 16 & 4 & 1\\\\1 & 161051 & 14641 & 1331 & 121 & 11 & 1\\\\1 & 1048576 & 65536 & 4096 & 256 & 16 & 1\\\\1 & 1889568 & 104976 & 5832 & 324 & 18 & 1\\\\1 & 3200000 & 160000 & 8000 & 400 & 20 & 1\\\\1 & 5153632 & 234256 & 10648 & 484 & 22 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1,    1024,    256,    64,  16,  4, 1],\n",
       "[1,  161051,  14641,  1331, 121, 11, 1],\n",
       "[1, 1048576,  65536,  4096, 256, 16, 1],\n",
       "[1, 1889568, 104976,  5832, 324, 18, 1],\n",
       "[1, 3200000, 160000,  8000, 400, 20, 1],\n",
       "[1, 5153632, 234256, 10648, 484, 22, 1]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = Matrix(np.array(X.speed)[:, np.newaxis]**[0, 5, 4, 3, 2, 1, 0])\n",
    "M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 37
    },
    "id": "BHktXffyzeYT",
    "outputId": "f1087528-a217-4eae-f521-b4addb204029"
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = M.pinv() * Y\n",
    "Y_ = M * b\n",
    "(Y - Y_).norm()"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "7.4.3項についての補足",
   "provenance": [],
   "toc_visible": true
  },
  "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
}