ex11notebookA.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "toc_visible": true,
      "authorship_tag": "ABX9TyNFEMNxwCGSebzrtB9MeZqm",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/takatakamanbou/MVA/blob/2023/ex11notebookA.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# MVA2023 ex11notebookA\n",
        "\n",
        "<img width=64 src=\"https://www-tlab.math.ryukoku.ac.jp/~takataka/course/MVA/MVA-logo11.png\"> https://www-tlab.math.ryukoku.ac.jp/wiki/?MVA/2023"
      ],
      "metadata": {
        "id": "P4preLw7SzqH"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# いつものいろいろインポート\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import pandas as pd\n",
        "import seaborn\n",
        "seaborn.set()\n",
        "\n",
        "# SciPy の 統計関数群の中の正規分布モジュール (scipy.stats.norm) と多変量正規分布モジュール (scipy.stats.multivariate_normal)\n",
        "from scipy.stats import norm, multivariate_normal\n",
        "\n",
        "# scikit-learn の線形判別分析と二次判別分析のクラス\n",
        "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis, QuadraticDiscriminantAnalysis\n",
        "\n",
        "# scikit-learn の主成分分析のクラス\n",
        "from sklearn.decomposition import PCA"
      ],
      "metadata": {
        "id": "9HCsHgFjR6jo"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# 人間 vs ほげ星人 (2)\n",
        "URL = 'https://www-tlab.math.ryukoku.ac.jp/~takataka/course/MVA/humanvshoge2.csv'\n",
        "dfHoge2 = pd.read_csv(URL)\n",
        "#dfHoge2"
      ],
      "metadata": {
        "id": "m6t8dQDbouxw"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# フィッシャーのアヤメのデータ\n",
        "from sklearn.datasets import load_iris\n",
        "iris = load_iris(as_frame=True)\n",
        "dfIris = iris.frame\n",
        "#dfIris"
      ],
      "metadata": {
        "id": "WlOIB6N3pFD8"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "## 判別分析 (2)\n",
        "---\n",
        "\n"
      ],
      "metadata": {
        "id": "Snu_5I6aGJS8"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "前回に続き，判別分析について解説します．"
      ],
      "metadata": {
        "id": "5kZCI7A5Q_Pg"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "### 多クラスの場合の判別分析\n",
        "\n",
        "前回は2クラスの問題に限定して解説していましたが，判別分析は多クラスの（クラス数3以上の）問題にも適用できます．\n"
      ],
      "metadata": {
        "id": "9xEU7w1Ic1NQ"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### 考え方"
      ],
      "metadata": {
        "id": "4t9_dz9dcN8R"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "データを「クラス$1$」から「クラス$K$」までの $K$ 通りのクラスに判別する問題の場合，次のようにします：\n",
        "1. クラス $k$ ($k = 1, 2, \\ldots, K$) のデータが正規分布 ${\\cal N}\\left(\\mathbf{\\mu}_k, \\Sigma_k\\right)$ に従うと仮定して，パラメータ $\\mathbf{\\mu}_k, \\Sigma_k$ を推定する．\n",
        "1. クラスが未知のデータ $\\mathbf{x}$ が与えられたら，それが ${\\cal N}\\left(\\mathbf{\\mu}_k, \\Sigma_k\\right)$ から得られたとすることの対数尤度 $\\log\\ell_k(\\mathbf{x})$ ($k = 1, 2, \\ldots, K$) を求め，どのクラスに対して尤度が最大となるか調べる：\n",
        "\n",
        "$$\n",
        "\\DeclareMathOperator*{\\argmax}{argmax}\n",
        "k^{*} = \\argmax_{k=1,\\ldots,K}{\\log\\ell_k(\\mathbf{x})} \\qquad (1)\n",
        "$$\n",
        "\n",
        "ここで，$\\displaystyle \\DeclareMathOperator*{\\argmax}{argmax} \\argmax_x f(x)$ は，$f(x)$ が最大となるときの $x$ の値を返す演算子です．$\\displaystyle \\DeclareMathOperator*{\\argmin}{argmin} \\argmin$ というのもあり，\n",
        "\n",
        "$$\n",
        "\\DeclareMathOperator*{\\argmax}{argmin}\n",
        "\\begin{aligned}\n",
        "\\min_{x}\\{(x-1)^2 + 2\\} &= 2\\\\\n",
        "\\argmin_{x}\\{(x-1)^2 + 2\\} &= 1\\\\\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "ということになります．つまり，式$(1)$ の $k^{*}$ は，尤度最大となったクラスの番号となります．"
      ],
      "metadata": {
        "id": "ehjr9ZTg9JDx"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "上記では，分散共分散行列がクラス毎に異なってもよいという設定で説明していますが，前回も説明したように，線形判別分析の場合は，全てのクラスの分散共分散行列が共通である，つまり，\n",
        "$\n",
        "\\Sigma_1 = \\Sigma_2 = \\dots = \\Sigma_K = \\Sigma\n",
        "$\n",
        "として考えることになります．"
      ],
      "metadata": {
        "id": "KDq6NfaTbJKA"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### 実験: Fisher のアヤメのデータの線形判別分析"
      ],
      "metadata": {
        "id": "pU2FDWNpcQ8G"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Fisher のアヤメのデータにはクラスが3つあります．線形判別分析してみましょう．"
      ],
      "metadata": {
        "id": "3SkQ2EsacYEa"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "## データの準備\n",
        "\n",
        "# 被説明変数（'target'列）を除いたものを np.array に\n",
        "X_iris = dfIris.drop(columns='target').to_numpy()\n",
        "print('# X_iris の最初の5行と X_iris.shape')\n",
        "print(X_iris[:5, :], X_iris.shape)\n",
        "print()\n",
        "\n",
        "# 'target' 列にクラス番号が格納されている\n",
        "Y_iris = dfIris['target'].to_numpy(dtype=int)\n",
        "for k, tn in enumerate(iris.target_names):\n",
        "    print(f'Class{k}: {tn}')\n",
        "print()\n",
        "print('# Y_iris と Y_iris.shape')\n",
        "print(Y_iris, Y_iris.shape)"
      ],
      "metadata": {
        "id": "4b_nzuX5c4FB"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# 線形判別分析\n",
        "lda = LinearDiscriminantAnalysis()\n",
        "lda.fit(X_iris, Y_iris)  # パラメータの推定\n",
        "Yt = lda.predict(X_iris) # 予測\n",
        "print('# クラス予測結果')\n",
        "print(Yt)\n",
        "ncorrect = np.sum(Yt == Y_iris)\n",
        "print(f'正解数/データ数 = {ncorrect}/{len(X_iris)}')"
      ],
      "metadata": {
        "id": "sdphH83eeemn"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "### 二次判別分析"
      ],
      "metadata": {
        "id": "pK6EK0RddMBL"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### 二次判別とは\n",
        "\n",
        "線形判別分析では，各クラスの正規分布の分散共分散行列が共通である，すなわち\n",
        "$\n",
        "\\Sigma_1 = \\Sigma_2 = \\dots = \\Sigma_K = \\Sigma\n",
        "$ が成り立つという仮定をおいていました．この仮定をなくして，分散共分散行列はクラスごとに異なってもよいとした場合を考えてみましょう．\n",
        "\n",
        "データを「クラス$1$」から「クラス$K$」までの $K$ 通りに判別する場合，データの所属するクラスを予測する手続きは次のようになります．\n",
        "\n",
        "1. クラス $k$ ($k = 1, 2, \\ldots, K$) のデータが正規分布 ${\\cal N}\\left(\\mathbf{\\mu}_k, \\Sigma_k\\right)$ に従うと仮定して，パラメータ $\\mathbf{\\mu}_k, \\Sigma_k$ を推定する．\n",
        "1. クラスが未知のデータ $\\mathbf{x}$ が与えられたら，それが ${\\cal N}\\left(\\mathbf{\\mu}_k, \\Sigma_k\\right)$ から得られたとすることの対数尤度 $\\log\\ell_k(\\mathbf{x})$ ($k = 1, 2, \\ldots, K$) を求め，尤度が最大となるクラス $k^{*}$ を見つける：\n",
        "\n",
        "$$\n",
        "\\log{\\ell_k(\\mathbf{x})} = -\\frac{D}{2}\\log(2\\pi) - \\frac{1}{2}\\log{|\\Sigma_k|} -\\frac{1}{2} (\\mathbf{x}-\\mathbf{\\mu}_k)^{\\top}\\Sigma_k^{-1}(\\mathbf{x}-\\mathbf{\\mu}_k)  \\qquad (k = 1, 2, \\ldots, K)\n",
        "$$\n",
        "\n",
        "$$\n",
        "\\DeclareMathOperator*{\\argmax}{argmax}\n",
        "k^{*} = \\argmax_{k=1,\\ldots,K}{\\log\\ell_k(\\mathbf{x})}\n",
        "$$\n",
        "\n",
        "$K = 2$ の場合，判別関数を $\\log\\ell_{1}(\\mathbf{x}) - \\log\\ell_{2}(\\mathbf{x})$ とすれば，その符号によってクラスの予測値を得ることができます（0以上ならクラス1，さもなくばクラス2）．ただし，判別関数が $\\mathbf{x}$ の一次式であった線形判別分析のときと異なり，この場合の判別関数は，$\\mathbf{x}$ の二次式となります．そのため，このように分散共分散行列がクラスごとに異なってもよいとした場合の判別分析のことを，**二次判別分析**(Quadratic Discriminant Analysis, QDA) といいます．\n",
        "\n",
        "判別の境界（判別関数の値が $0$ となる点の集合）は，線形判別では平面でしたが，二次判別では二次曲面となります．\n"
      ],
      "metadata": {
        "id": "Aq6UFH0f0ygb"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### 例: 「人間」vs「ほげ星人」の二次判別\n",
        "\n",
        "2次元2クラスのデータを使って実際にやってみましょう．"
      ],
      "metadata": {
        "id": "yjGRYkS94vyv"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# データの準備\n",
        "X_hoge = dfHoge2[['Height', 'Weight']].to_numpy()\n",
        "Y_hoge = (dfHoge2['Class'] == 'Human').to_numpy(dtype=int)\n",
        "print(X_hoge[:5, :], X_hoge.shape)\n",
        "print(Y_hoge, Y_hoge.shape)"
      ],
      "metadata": {
        "id": "MMyb_MxRxnD7"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# 散布図\n",
        "xmin, xmax = 0, 250\n",
        "ymin, ymax = 0, 150\n",
        "fig, ax = plt.subplots(figsize=(12, 6))\n",
        "ax.scatter(X_hoge[Y_hoge == 1, 0], X_hoge[Y_hoge == 1, 1], label='Human')\n",
        "ax.scatter(X_hoge[Y_hoge == 0, 0], X_hoge[Y_hoge == 0, 1], label='Hoge')\n",
        "ax.set_xlim(xmin, xmax)\n",
        "ax.set_xlabel('Height')\n",
        "ax.set_ylim(ymin, ymax)\n",
        "ax.set_ylabel('Weight')\n",
        "ax.set_aspect('equal')\n",
        "ax.legend()\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "vjQEVO4dxhCp"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "このデータは，前回使っていた「人間」vs「ほげ星人」のデータとは別物です．\n",
        "前回使ったものでは2つのクラスの分布が似た形をしていましたが，こちらでは，両者の分布の形が少し違っています．線形判別と二次判別の結果の違いを示すため，こちらのデータを使ってみることにします（前回のデータだと，線形判別でも二次判別でもほとんど同じ結果となります）．"
      ],
      "metadata": {
        "id": "n_cDADg5yJnw"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# 線形判別分析モデルのパラメータの推定\n",
        "lda = LinearDiscriminantAnalysis(store_covariance=True)\n",
        "lda.fit(X_hoge, Y_hoge)\n",
        "\n",
        "# 二次判別分析モデルのパラメータの推定\n",
        "qda = QuadraticDiscriminantAnalysis(store_covariance=True)\n",
        "qda.fit(X_hoge, Y_hoge)"
      ],
      "metadata": {
        "id": "AX6uKq84y70h"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "fig, ax = plt.subplots(1, 2, figsize=(16, 8))\n",
        "\n",
        "## 確率密度と判別結果の描画のためのグリッドデータの準備\n",
        "xmin, xmax = 0, 250\n",
        "ymin, ymax = 0, 150\n",
        "xx, yy = np.mgrid[xmin:xmax:1, ymin:ymax:1]\n",
        "XX = np.dstack((xx, yy))\n",
        "\n",
        "##  散布図\n",
        "for i in [0, 1]:\n",
        "    ax[i].scatter(X_hoge[Y_hoge == 1, 0], X_hoge[Y_hoge == 1, 1], label='Human')\n",
        "    ax[i].scatter(X_hoge[Y_hoge == 0, 0], X_hoge[Y_hoge == 0, 1], label='Hoge')\n",
        "    ax[i].set_xlim(xmin, xmax)\n",
        "    ax[i].set_xlabel('Height')\n",
        "    ax[i].set_ylim(ymin, ymax)\n",
        "    ax[i].set_ylabel('Weight')\n",
        "    ax[i].set_aspect('equal')\n",
        "    ax[i].legend()\n",
        "\n",
        "## 確率密度の等高線\n",
        "# LDA\n",
        "zz = multivariate_normal.pdf(XX, mean=lda.means_[1], cov=lda.covariance_)\n",
        "ax[0].contour(xx, yy, zz, colors='blue')\n",
        "zz = multivariate_normal.pdf(XX, mean=lda.means_[0], cov=lda.covariance_)\n",
        "ax[0].contour(xx, yy, zz, colors='red')\n",
        "# QDA\n",
        "zz = multivariate_normal.pdf(XX, mean=qda.means_[1], cov=qda.covariance_[1])\n",
        "ax[1].contour(xx, yy, zz, colors='blue')\n",
        "zz = multivariate_normal.pdf(XX, mean=qda.means_[0], cov=qda.covariance_[0])\n",
        "ax[1].contour(xx, yy, zz, colors='red')\n",
        "\n",
        "## 判別結果で領域を塗り分ける\n",
        "# LDA\n",
        "zz = lda.predict(XX.reshape(-1, 2)).reshape((XX.shape[0], XX.shape[1]))\n",
        "ax[0].contourf(xx, yy, zz,   cmap='Blues',   alpha=0.2)\n",
        "ax[0].contourf(xx, yy, 1-zz, cmap='Oranges', alpha=0.2)\n",
        "# QDA\n",
        "zz = qda.predict(XX.reshape(-1, 2)).reshape((XX.shape[0], XX.shape[1]))\n",
        "ax[1].contourf(xx, yy, zz,   cmap='Blues',   alpha=0.2)\n",
        "ax[1].contourf(xx, yy, 1-zz, cmap='Oranges', alpha=0.2)\n",
        "\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "Sr9Xl4HQzMEN"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "左が線形判別，右が二次判別の結果です．線形判別では，2クラスの分散共分散行列が共通であると仮定していますので，描かれた2つの同心楕円は同じ形をしています．一方，二次判別では両者は別々ですので，異なる形となっています．\n",
        "\n",
        "この図には，判別の境界も描いてあります．線形判別では直線であるのに対して，二次判別では放物線となっています．\n",
        "二次判別で得られる2クラスの境界は，この例では放物線ですが，放物線（面）以外の二次曲線（面）になることもあります（例: https://scikit-learn.org/stable/modules/lda_qda.html ）．\n"
      ],
      "metadata": {
        "id": "NqhcccTu1jfF"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "### 補足"
      ],
      "metadata": {
        "id": "WsqwZBErGt5h"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### データの標準化"
      ],
      "metadata": {
        "id": "YlvlpHD7n1AP"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "「主成分分析(2)」の回の「主成分分析とデータの標準化」という節で，データの標準化が主成分分析の結果に与える影響について考えました．\n",
        "データによっては変数間で値の範囲・大きさ範囲が全く異なることがあるので，データを標準化してから分析した方がよい場合がある，ということでした．実例は示しませんが，判別分析においても同様のことが言えます．"
      ],
      "metadata": {
        "id": "mwLtBwi3qrDY"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### 線形判別 vs 二次判別／正規分布の仮定が成り立たないと？\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "1PzhALiQ5x06"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "線形判別分析では判別の境界が平面でしたが，二次判別分析では二次曲面になります．そのため，クラス毎の分散共分散行列が共通であるという条件が成り立たないような場合には，二次判別の方が良い（予測の精度が高い）結果が得られる可能性があります．しかし，実際には，次のようなことに注意しなければなりません．\n",
        "\n",
        "- 二次判別の方が推定すべきパラメータの数が多い： 分散共分散行列ひとつは $\\frac{1}{2}D(D+1)$ 個のパラメータを持ちます．二次判別の場合，これをクラスごとに推定しますので，パラメータ数がさらに $K$ 倍になります．サンプルサイズが小さいデータだと，分散共分散行列の推定精度が悪くて良い結果が得られないことがあります．\n",
        "- 正規分布の仮定が成り立たない場合，線形判別でも二次判別でも良い結果は得られないことがある： 線形判別も二次判別も，クラスごとのデータが正規分布に従うことを仮定しています．この仮定が満たされないような場合には，どちらの手法でも良い結果が得られないことがあります．\n",
        "\n",
        "次の例は，データの分布が明らかに正規分布ではないような場合にどうなるかを示しています．"
      ],
      "metadata": {
        "id": "zautJh1XDVLe"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "### 正規分布に従わない人工データ\n",
        "from sklearn.datasets import make_moons\n",
        "X_moon, Y_moon = make_moons(n_samples=300, noise=0.2)\n",
        "\n",
        "### QDA してみる\n",
        "qda = QuadraticDiscriminantAnalysis(store_covariance=True)\n",
        "qda.fit(X_moon, Y_moon)\n",
        "\n",
        "### 散布図と判別結果\n",
        "fig, ax = plt.subplots(1, 2, figsize=(16, 8))\n",
        "\n",
        "## 確率密度と判別結果の描画のためのグリッドデータの準備\n",
        "xmin, xmax = -1.5, 2.5\n",
        "ymin, ymax = -1.0, 1.5\n",
        "xx, yy = np.mgrid[xmin:xmax:0.01, ymin:ymax:0.011]\n",
        "XX = np.dstack((xx, yy))\n",
        "\n",
        "##  散布図\n",
        "for i in [0, 1]:\n",
        "    ax[i].scatter(X_moon[Y_moon == 1, 0], X_moon[Y_moon == 1, 1], label='Class1')\n",
        "    ax[i].scatter(X_moon[Y_moon == 0, 0], X_moon[Y_moon == 0, 1], label='Class0')\n",
        "    ax[i].set_xlim(xmin, xmax)\n",
        "    ax[i].set_ylim(ymin, ymax)\n",
        "    ax[i].set_aspect('equal')\n",
        "    ax[i].legend()\n",
        "\n",
        "## 確率密度の等高線\n",
        "zz = multivariate_normal.pdf(XX, mean=qda.means_[1], cov=qda.covariance_[1])\n",
        "ax[1].contour(xx, yy, zz, colors='blue')\n",
        "zz = multivariate_normal.pdf(XX, mean=qda.means_[0], cov=qda.covariance_[0])\n",
        "ax[1].contour(xx, yy, zz, colors='red')\n",
        "\n",
        "## 判別結果で領域を塗り分ける\n",
        "zz = qda.predict(XX.reshape(-1, 2)).reshape((XX.shape[0], XX.shape[1]))\n",
        "ax[1].contourf(xx, yy, zz,   cmap='Blues',   alpha=0.2)\n",
        "ax[1].contourf(xx, yy, 1-zz, cmap='Oranges', alpha=0.2)\n",
        "\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "t1ZzAl5C2kai"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "当たり前ですが，あまりうまく判別できていません．\n",
        "\n",
        "「機械学習I/II」では，このような場合でもうまく判別できるような手法を取り上げます．"
      ],
      "metadata": {
        "id": "yDOdYtGKBgWN"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### 事前確率（クラスごとのデータの出現確率）の扱い"
      ],
      "metadata": {
        "id": "amEx2o0eDiWd"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "ここまで解説してきた判別分析の方法は，クラスごとの正規分布を推定し，それらの対数尤度を基準として判別するというものです．「クラスごとのデータが正規分布に従う」という仮定をおいており，その仮定が成り立たないときにはうまくいかないかもしれない，ということはすでに説明しました．\n",
        "\n",
        "実は，ここまで説明してきた方法には，隠れた仮定がもうひとつあります．それは，「どのクラスのデータも等しい確率で生起する」というものです．この仮定が成り立たない場合，すなわち，クラスごとのデータの出現確率に偏りがあるような場合，クラスごとの正規分布に対する対数尤度を基準とする，という方法ではよい結果が得られない可能性があります．\n",
        "\n",
        "これに対処する方法の一つが，クラスごとのデータの出現確率（これを「クラス事前確率」といいます）まで含めてモデル化を行う，というものです．\n",
        "具体的な定式化等については，「機械学習I/II」で取り上げます．scikit-learn の判別分析のクラスでは，このクラス事前確率を含めたモデルが使われており，データ中の各クラスの出現頻度からその値を推定するようになっています（注）．\n",
        "\n",
        "※注意: 以下のリンク先のドキュメントの $P(y=k)$ というのが事前確率です：\n",
        "https://scikit-learn.org/stable/modules/lda_qda.html#mathematical-formulation-of-the-lda-and-qda-classifiers"
      ],
      "metadata": {
        "id": "V4lnJjl4DiWl"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "### 【発展】次元削減手法としての線形判別分析"
      ],
      "metadata": {
        "id": "wmLgY-V0dCAs"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "線形判別分析は判別のための手法ですが，実は，次元削減手法の一種として考えることもできます．ここでは，詳しい解説は省いて，直感的にそのことを説明します．"
      ],
      "metadata": {
        "id": "Mz_Y-qoXhbUM"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### 2クラス線形判別はデータを1次元に変換する"
      ],
      "metadata": {
        "id": "G65F5RgWeiVF"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "2クラス問題に対する線形判別分析は，データを1次元に次元削減する軸を見つける処理とみなせます．\n",
        "\n",
        "クラス1に属するデータが ${\\cal N}(\\mathbf{\\mu}_1, \\Sigma)$ に従い，クラス2に属するデータが ${\\cal N}(\\mathbf{\\mu}_2, \\Sigma)$ に従うとき，線形判別分析によって得られる線形判別関数 $z(\\mathbf{x})$ は\n",
        "\n",
        "$$\n",
        "z(\\mathbf{x}) = (\\mathbf{\\mu}_1 - \\mathbf{\\mu}_2)^{\\top}\\Sigma^{-1}\\left( \\mathbf{x} - \\frac{\\mathbf{\\mu}_1+\\mathbf{\\mu}_2}{2} \\right)\n",
        "$$\n",
        "\n",
        "と表せるのでした（ex10notebookAの式$(16)$）．\n",
        "ここで $\\mathbf{w} = \\Sigma^{-1}(\\mathbf{\\mu}_1 - \\mathbf{\\mu}_2)$ とおくと，\n",
        "\n",
        "$$\n",
        "z(\\mathbf{x}) = \\mathbf{w}^{\\top}\\left( \\mathbf{x} - \\frac{\\mathbf{\\mu}_1+\\mathbf{\\mu}_2}{2} \\right)\n",
        "$$\n",
        "\n",
        "と書けます．この式は，$D$ 個の要素から成る $\\mathbf{x}$ から一つの実数値 $z(\\mathbf{x})$ を求める変換を表しています（注）．$z(\\mathbf{x})$ の値は，下図左に示すように，点 $\\frac{\\mathbf{\\mu}_1+\\mathbf{\\mu}_2}{2}$ を通りベクトル $\\mathbf{w}$ に並行な直線上に点 $\\mathbf{x}$ を射影したときに，この直線上で測ったその点の位置を表しています．下図右は，$z(\\mathbf{x})$ の値を横軸にとってヒストグラムを描いたものです．この軸上で見ると，点Aは正の方に，点Bは負の方にいることになります．\n",
        "\n",
        "※注意: 主成分分析において $\\mathbf{x}$ から主成分スコア $y$ を求める計算（$y = \\mathbf{u}^{\\top}(\\mathbf{x}-\\bar{\\mathbf{x}}）$) と似た形をしてますね．"
      ],
      "metadata": {
        "id": "cyIHOqsfmE9Q"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<img src=\"https://www-tlab.math.ryukoku.ac.jp/~takataka/course/MVA/dimreductionLDA.png\">"
      ],
      "metadata": {
        "id": "zlZHsnM4xkQn"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "2クラス問題の線形判別分析でデータを1次元に変換できることが分かりましたが，この変換は，主成分分析で得られるものとはどう違うのでしょうか．両者の違いを理解してもらうために，2クラス2次元のデータに対して主成分分析と線形判別分析を適用した結果を例にあげます．\n",
        "\n",
        "下図左のピンク色の星印は，これらデータ全体の（クラスを区別しないで求めた）平均です．このデータでは，2クラスに等しい数のデータが含まれているので，これは2クラス（オレンジ色と青色）それぞれの平均の中点と一致しています．\n",
        "\n",
        "このデータに主成分分析を適用して第1主成分を求めると，ピンクの矢線の向きとなります．主成分分析ではデータのクラス分けは考慮されず，データ全体の分散が最も大きくなる軸が選ばれています．一方，\n",
        "線形判別分析を適用して得られる軸は水色の矢線の向きとなります．\n",
        "\n",
        "どちらの軸を選んでも2次元のデータを1次元にすることができますが，その軸方向に見たデータの分布は異なります．主成分分析を用いて1次元に変換したデータの値（第1主成分スコア）は，下図真ん中のように分布しており，2つのクラスの値が結構入り混じっています．\n",
        "これに対して，線形判別分析を用いて1次元に変換したデータの値（判別関数の値）は，下図右のように分布します．こちらの方が2クラスのデータがよりうまく分離されていることが分かります．\n"
      ],
      "metadata": {
        "id": "GlH4_KO_ef0s"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<img src=\"https://www-tlab.math.ryukoku.ac.jp/~takataka/course/MVA/PCAvsLDA.png\">"
      ],
      "metadata": {
        "id": "v4T_hlMEd16P"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "このように，線形判別分析を利用した次元削減では，クラスをよく分離するような軸が選ばれ，主成分分析を利用した次元削減では，データ全体の分散が大きくなるような軸が選ばれます．判別が目的である場合，主成分分析による次元削減は適切ではないこともあります．\n"
      ],
      "metadata": {
        "id": "Bpisb2K8pNfo"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### $K$クラス線形判別による $(K-1)$ 次元への次元削減"
      ],
      "metadata": {
        "id": "9awDmcPUmlYB"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "上記の説明は，2クラス問題のデータを1次元に変換するという話でした．これを一般化すると，次のことがいえます．\n",
        "\n",
        "> $K$クラス問題に線形判別分析を適用すると，データを $(K-1)$ 次元に変換する線形変換の中で，クラス間の分離が最もよくなるようなものを見つけることができる．\n",
        "\n",
        "このことは，数学的にきちんと証明できるのですが，ここではその証明およびこの線形変換の具体的な計算法についての説明は省略します（注）．\n",
        "変換後の次元数を元の次元数以下で自由に決められる主成分分析と違い，次元数が $(K-1)$ 以下になる（$K$次元目以降は有効ではない）という強い制約があることに注意が必要です．\n",
        "\n",
        "<span style=\"font-size: 75%\">\n",
        "※注: 気になるひとは，この授業の「参考情報」のページに上げた「わかりやすいパターン認識 第2版」をどうぞ．\n",
        "</span>"
      ],
      "metadata": {
        "id": "tfPkhgyhpS2_"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "試しに，Fisherのアヤメのデータを次元削減する実験をやってみましょう．このデータは下図に示すような4次元のデータです．クラス数は3です．主成分分析と線形判別分析のそれぞれを利用して，このデータを2次元に変換してみましょう．"
      ],
      "metadata": {
        "id": "b4kuc7Phpukq"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# 4つの変数のうち2つ選んで散布図を描く\n",
        "seaborn.pairplot(dfIris, hue='target', palette='tab10', corner=True)"
      ],
      "metadata": {
        "id": "lN85T-nTpf8c"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "まずは主成分分析です．以前の notebook でやっていたように NumPy 等を使って自分で計算してもいいのですが，ここでは scikit-learn の [sklearn.decomposition.PCA](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html) を使っています．主成分分析ではクラス分けを考慮しませんから，`X` のみを用い，`Y` の方は使いません．\n",
        "\n"
      ],
      "metadata": {
        "id": "eusfdrlOsFfq"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# 主成分分析による次元削減\n",
        "pca = PCA(n_components=2)\n",
        "pca.fit(X_iris)\n",
        "X_pca = pca.transform(X_iris)\n",
        "print(X_pca[:5, :], X_pca.shape)"
      ],
      "metadata": {
        "id": "p9mLWwPdrCQb"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "`X_pca` の0列目が第1主成分スコア，1列目が第2主成分スコア，です．\n",
        "\n",
        "次に，線形判別分析による次元削減です．前回の演習でも使っていた [sklearn.discriminant_analysis.LinearDiscriminantAnalysis](https://scikit-learn.org/stable/modules/generated/sklearn.discriminant_analysis.LinearDiscriminantAnalysis.html) に次元削減の変換を行う機能も備わっています．"
      ],
      "metadata": {
        "id": "od6QF2LKs3Nf"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "左が主成分分析を利用した変換（上位2つの主成分を用いて2次元に次元削減したもの）の結果，右が線形判別分析を利用した変換の結果です．主成分分析で2次元にした場合，3つのクラスがあまりよく分離されていません．しかし，線形判別分析で2次元にした場合，3つのクラスがよく分離されています（横軸の値だけでほぼ3クラスが判別できる）．"
      ],
      "metadata": {
        "id": "Me5NY6-7uTl7"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# 線形判別分析による次元削減\n",
        "lda = LinearDiscriminantAnalysis()\n",
        "lda.fit(X_iris, Y_iris)\n",
        "X_lda = lda.transform(X_iris)\n",
        "print(X_lda[:5, ], X_lda.shape)"
      ],
      "metadata": {
        "id": "rgqIcwXLqkqO"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "元のアヤメのデータは4次元ですが，2次元に変換すれば散布図に描いてみることができます．描いてみるとこんなんなります．"
      ],
      "metadata": {
        "id": "cs2jTBqqt9tD"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "fig, ax = plt.subplots(1, 2, figsize=(16, 8))\n",
        "\n",
        "for k in range(3):\n",
        "    ax[0].scatter(X_pca[Y_iris == k, 0], X_pca[Y_iris == k, 1], label=iris.target_names[k])\n",
        "    ax[1].scatter(X_lda[Y_iris == k, 0], X_lda[Y_iris == k, 1], label=iris.target_names[k])\n",
        "\n",
        "for i in [0, 1]:\n",
        "    ax[i].axhline(0, color='gray')\n",
        "    ax[i].axvline(0, color='gray')\n",
        "    ax[i].set_aspect('equal')\n",
        "    ax[i].legend()\n",
        "\n",
        "ax[0].set_xlim(-4.5, 4.5)\n",
        "ax[0].set_ylim(-2.5, 2.5)\n",
        "ax[1].set_xlim(-10.5, 10.5)\n",
        "ax[1].set_ylim(-6, 6)\n",
        "\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "dyUwSp07rSEN"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "Ij6JZYwaE8BA"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}