「入門はじめての多変量解析」をPythonで写経 Vol.22 ～ ５章「はじめての判別分析」①線形判別分析・プロローグ

５章「はじめての判別分析」

書籍の著者　石村貞夫 先生、石村光資郎 先生

書籍「入門はじめての多変量解析」５章「はじめての判別分析」の Python写経活動記録 です。 

多変量解析の入門を Python と一緒に学ぶ写経シリーズです。

判別分析はいわゆる「分類」タスクです。
書籍は次の２種類の判別分析を掲載しています。
1️⃣ 線形判別分析（２変数、２群） ◀ 今回はこちら！
2️⃣ マハラノビス距離による判別

この記事は、線形判別分析をざっくり実践 して、判別分析の外観を体感します。
Python の２つのライブラリで線形判別分析を試してみます！

ChatGPT 活用型学習で進めてまいります！
では書籍を開いて多変量解析の旅に出かけましょう🚀

はじめに

このブログシリーズは、書籍「入門はじめての多変量解析」（東京図書、「テキスト」と呼びます）の Python 写経を通じて得た「多変量解析の楽しさ」をご紹介します。

書籍の紹介と引用表記はリンク先の記事に掲載しています。

５章 はじめての判別分析

この記事は５章の以下のSectionを取り扱います。

記事に用いるデータは、テキストに掲載されたデータそのものを引用しています。
データ件数の少ないものはコード上でデータを登録し、データ件数の多いものはCSVファイル化してデータを読み込みしています。

この記事で用いるライブラリをインポートします。

### インポート

# 数値計算
import numpy as np
import pandas as pd

# 線形判別分析
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from mlxtend.feature_extraction import LinearDiscriminantAnalysis as lda_mlxtend

# 決定境界描画
from sklearn.inspection import DecisionBoundaryDisplay
from matplotlib import colors

# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

イントロダクション

判別分析の超イントロをChatGPTに書いてもらいました。
どうぞご覧ください！

想像してみてください。あなたの手元に、いくつかの果物が並んでいるとします。リンゴとオレンジをぱっと見分けるのは簡単ですよね😊
では、その「見た目」「重さ」「色」の情報をコンピュータに教えて、初めて見る果物がリンゴかオレンジか自動で判別してくれたら面白くないですか？

そんな魔法のような仕組みが、統計の「判別分析」です。

• 何をする？
  過去のデータ（例：色・形・重さ）をもとに、ものごとがどのグループに属するかを見分けるルールを学び取ります。

• どんなとき？
  新しいお客さんが「どの会員ランクに当てはまるか」知りたいとき。
  メールが「迷惑メールかどうか」を振り分けたいとき。

• どうして楽しい？
  自分のデータに「仕分け機」をつくる感覚！
  一度ルールを覚えさせれば、あとは自動でジャッジしてくれるので、まるで自分専用の名札発行マシンができたようなワクワクが味わえます✨

この記事では、判別分析の基本イメージからはじめて、実際にデータを使って「どこで線を引くか」を学んでいきます。
手を動かすだけで、データの世界が一気に身近になりますよ！

さあ、一緒に判別の魔法を体験してみましょう🔍🎉

データを仕分けする、名札を付ける…　いい表現ですね！

テキストでは主に次のことを学べます。
・２つのグループに仕分けする「二値分類」
・グループを仕分けする「境界線」（直線／曲線）を見つける

では、テキストの例題データをお借りして、判別分析の世界を「Python」でざっくり体感しに行きましょう！

判別分析の対象データ p.178～

テキスト p.179 表 5.1.1「前立腺疾患のデータ」をお借りします。
グループ１とグループ２で構成されています。
このグループ分けが「判別・分類」の対象になります。

こちらはグループ１（前立腺がん）です。

### マーカー測定結果 p.179 表5.1.1

# 前立腺ガンのグループG1
group1 = pd.DataFrame(
    {'被験者No.': range(1, 8),
     'マーカーA': [3.4, 3.9, 2.2, 3.5, 4.1, 3.7, 2.8],
     'マーカーB': [2.9, 2.4, 3.8, 4.8, 3.2, 4.1, 4.2],
     '結果': np.ones(7).astype(int)})
group1

【実行結果】
７名の測定データです。

続いてグループ２（前立腺肥大症）です。

# 前立腺肥大症のグループG2 p.179 表5.1.1
group2 = pd.DataFrame(
    {'被験者No.': range(8, 16),
     'マーカーA': [1.4, 2.4, 2.8, 1.7, 2.3, 1.9, 2.7, 1.3],
     'マーカーB': [3.5, 2.6, 2.3, 2.6, 1.6, 2.1, 3.5, 1.9],
     '結果': np.zeros(8).astype(int)})
group2

【実行結果】
８名の測定データです。

２つのグループを「結果」で区別しています。
・グループ１：結果は１（前立腺がん）
・グループ２：結果は０（前立腺肥大症）

２つのデータを結合します。

# データの結合
data1 = pd.concat([group1, group2], axis=0).reset_index(drop=True)
data1

【実行結果】

テキストは５章全体を通じて、マーカーＡ、マーカーＢの測定値からグループ１or２を判別するモデルを検討し、新しいデータを用いた予測に取り組みます。
検査結果から症状を予測するモデルを作る感じです。
そして次のデータが、予測に用いる新しいデータ（ S さんのデータ）です。

モデル構築しやすいようにデータセットを整形します。
上記データを「説明変数 X」と「目的変数 y」に分離します。
あわせて、Sさんデータを X_new で作成します。

### データセットの作成

# 学習データ
X = data1[['マーカーA', 'マーカーB']].values
y = data1['結果'].values

# 予測データ shape(n, d)
X_new = np.array([[2.7, 3.1]])

準備万端です！

線形判別関数 p.180～

可視化の力を借りて、判別・分類のイメージを膨らませていきましょう。

最初にデータの散布図を描きます。
グループごとに色を変えて可視化します。
seaborn ライブラリの可視化関数で 引数 hue を使うと、「結果」の値に応じた色分けができます。

### ２つのグループの散布図 p.180 図5.2.1

# 描画領域の設定
plt.figure(figsize=(6, 6))
# hue:結果に症名を設定
hue = ['グループG1' if data1['結果'].values[i]==1 else 'グループG2'
       for i in range(len(data1))]
# data1の散布図の描画
sns.scatterplot(data=data1, x='マーカーA', y='マーカーB', s=100,
                hue=hue, palette=['tomato', 'tab:blue'])
# Sさんの散布図の描画
sns.scatterplot(x=X_new[:, 0], y=X_new[:, 1], marker='*', s=500, color='orange',
                label='S')
# 修飾
plt.xlim(0, 5)
plt.ylim(1, 5)
plt.yticks(range(1, 6))
plt.gca().set_aspect('equal')
plt.legend(loc='upper left');

【実行結果】
グループ１が赤色、グループ２が青色、Sさんが星印です。

グループ１は、マーカーA、マーカーBの両方の値が大きい感じがします。
グループ２は、マーカーA、マーカーBの両方の値が小さい感じがします。

ええいっ！と大まかに区分けしてみましょう。

### ２つのグループの散布図 p.180 図5.2.1

# 描画領域の設定
plt.figure(figsize=(6, 6))
# hue:結果に症名を設定
hue = ['グループG1' if data1['結果'].values[i]==1 else 'グループG2'
       for i in range(len(data1))]
# data1の散布図の描画
sns.scatterplot(data=data1, x='マーカーA', y='マーカーB', s=100,
                hue=hue, palette=['tomato', 'tab:blue'])
x_vals = np.linspace(data1['マーカーA'].min(), data1['マーカーA'].max(), 101)
# Sさんの散布図の描画
sns.scatterplot(x=X_new[:, 0], y=X_new[:, 1], marker='*', s=500, color='orange',
                label='S')
# 境界線（仮）の描画
plt.plot(x_vals, 9.4 -2.4 * x_vals, color='tab:red', lw=1, ls='--')
# 修飾
plt.xlim(0, 5)
plt.ylim(1, 5)
plt.yticks(range(1, 6))
plt.gca().set_aspect('equal')
plt.legend(loc='upper left');

【実行結果】
適当な場所に境界線を引いて区切ってみました。
まあまあいい感じにグループ１とグループ２を仕分け＝判別できているように見えます。

上の図のように「直線的な境界線」を引くイメージで判別をするモデル式が 線形判別関数 $${z}$$ なのです。
テキストで扱う「２変数」の場合、線形判別関数は次の数式で示されます。

$${z=0}$$ となる $${0 =a_1 x_1 + a_2 x_2 + a_0 }$$ が、上のチャートで表現したい境界線の数式です。

線形判別関数を使うと…

という感じで数式計算で判別できるようになります！

線形判別関数の求め方はテキスト Section 5.3 ～ 5.5 でゆっくり・たっぷり解説されています。
次回記事で検討します！

この先は Python を用いて、判別分析のさまざまな出力結果を予習いたしましょう。
そして、線形判別の流派にはさまざまあり、テキストの線形判別関数は流派の１つであること（唯一の線形判別関数ではないこと）を体験します。

テキストの線形判別関数を再現するPythonコード

テキストの線形判別関数を再現できるように「イメージ」して、線形判別分析クラスを書きました。
今回の例題データに関しては、テキストでやりたいことをほぼ再現できています！
ただし、テキストと計算途中の小数点の丸め桁数が異なるため、テキストの計算結果と少し違います。

0️⃣ 準備
ひとまずオリジナルの線形判別分析クラス＆ヘルパー関数を定義します。
線形判別分析の実装の際、次のWebサイトの数式を参考にいたしました。
ありがとうございます！

## 線形判別分析クラスの定義
# 参考URL:
# https://hkawabata.github.io/technical-note/note/ML/Preprocess/lda.html

class LDA:

    # 初期化【引数】 n_discriminants: 判別子の数
    def __init__(self, n_discriminants=1):
        self.n_discriminants = n_discriminants

    # 学習関数【引数】X: 説明変数（numpy配列）、y:目的変数（numpy配列）
    def fit(self, X, y):

        ## 準備
        # 平坦な配列を列ベクトルに変換する関数の定義
        vec = lambda x: x.reshape(-1, 1)
        # 説明変数・目的変数をクラス属性に設定
        self.X, self.y = X, y

        ## 設定
        # クラスの要素, クラスの標本サイズ
        classes, ncs = np.unique(self.y, return_counts=True)
        self.classes = classes

        ## 平均の算出
        # 学習データの変数ごとの平均ベクトルm shape=(d, 1)
        self.means = self.X.mean(axis=0)
        m = vec(self.means)
        # クラスごとの各変数の平均行列 shape=(n_classes, d)
        mcs = np.vstack(
            [self.X[np.where(self.y == c)].mean(axis=0) for c in self.classes])

        ## SB, SW の算出、クラス属性に設定
        # グループ間変動SBの算出 shape=(d, d)
        # テキストSection5.4のグループ間変動に対応
        self.SB = sum(np.array(
            [nc * (vec(mc) - m) @ (vec(mc) - m).T for nc, mc in zip(ncs, mcs)]
        ))

        # グループ内変動SWの算出 shape=(d, d)
        # テキストSection5.4のグループ内変動に対応
        self.SW = sum(np.array(
            [(vec(xi) - vec(mc)) @ (vec(xi) - vec(mc)).T
             for c, mc in zip(self.classes, mcs)
             for xi in self.X[np.where(self.y == c)]]
        ))

        ## 固有値・固有ベクトルの算出
        e_vals, e_vecs = np.linalg.eig(np.linalg.inv(self.SW) @ self.SB)

        ## 固有値・固有ベクトルを固有値の昇順で並び替えて、クラス属性に設定
        sort_index = e_vals.argsort()[::-1]
        self.e_vals = e_vals[sort_index]
        self.e_vecs = e_vecs[:, sort_index]

    # 変換関数（固有ベクトルを用いる次元削減） ※mlxtendと一致する
    def transform(self, X):
        # データを固有ベクトルで変換
        X_tranformed = X @ self.e_vecs[:, :self.n_discriminants]
        # 戻り値：変換後のデータ
        return X_tranformed
    
    # 【変数が２つの場合】判別得点算出関数
    def discriminant_score_2d(self, X_new):
        # 係数の取得
        w = self.e_vecs[:, 0]
        # 傾きの算出
        slope = w[0] / w[1]
        # 切片の算出
        intercept = - (slope * self.means[0] + self.means[1])
        # 判別得点の算出：線形判別関数 z = a1*x1 + x2 + a0 で計算
        z = slope * X_new[:, 0] + X_new[:, 1] + intercept
        # 戻り値：判別得点z
        return z


## 【変数が２つの場合】ヘルパー関数の定義

# 【変数が２つの場合】判別得点算出関数（fit済みのldaインスタンスが必要）
def discriminant_score_2d(X_new, X_train, trans_matrix):
    # x2 = slope * x1 + interceptの傾きと切片を算出
    slope_, intercept_ = coef_2d(X_train, trans_matrix)
    # 傾きと切片の符号を逆転
    slope, intercept = -slope_, -intercept_
    # 判別得点の算出：線形判別関数 z = a1*x1 + x2 + a0 で計算
    z = slope * X_new[:, 0] + X_new[:, 1] + intercept
    return z

# 【変数が２つの場合】平面グラフに線形判別関数を描画する際の傾きと切片を求める関数
def coef_2d(X, trans_matrix):
    # 変数ごとの平均値の算出
    means = X.mean(axis=0)
    # 係数の取得
    w = trans_matrix[:, 0]
    # 傾きの算出
    slope = - w[0]/ w[1]
    # 切片の算出
    intercept = -slope * means[0] + means[1]
    # 戻り値：傾き、切片
    return slope, intercept

【実行結果】なし

使い方を確認しつつ、線形判別分析を進めていきましょう！

1️⃣ 線形判別分析の実行
某ライブラリっぽく、インスタンス生成とモデルの学習の２ステップで線形判別分析を実行します。

## 線形判別分析の実行

# 線形判別分析のインスタンス生成
lda = LDA(n_discriminants=1)
# モデルの学習
lda.fit(X, y)

【実行結果】なし
引数 n_discriminants=1 は、変換（次元削減）時の出力次元数です。
今回は１次元です。

2️⃣ 判別得点の算出
判別得点は 線形判別関数 $${z}$$ に変数 $${x_1, x_2}$$ の値を入れて計算した値です。
判別得点$${>0}$$ の場合はグループ１（前立腺がん）、$${<0}$$ の場合はグループ２（前立腺肥大症）を示します。

学習済みのモデルに対して「discriminant_score_2d」メソッドを適用します。
引数で判別得点を算出したいデータを与えます。
今回は学習データ X の判別得点を算出します。

## 判別得点の算出：メソッド版
lda_score = lda.discriminant_score_2d(X)
lda_score

【実行結果】
被験者No. １～15 の判別得点です。
テキスト p.202 表 5.6.1 と表 5.6.2 に相当します。

前半７人がグループ１でプラスになるはず、後半８人がグループ２でマイナスになるはず、です。
14 番目の判別得点がプラスの値になっていて、グループ１だと誤判別しています。

テキストの計算結果との違いは、計算途中の小数点の丸め桁数の違いによるものです。
結構なずれが発生するものなのですね。

◆ 線形判別関数の確認
線形判別関数の傾き、切片を表示しましょう。

# 傾き、切片の表示
print('傾き a1: ', lda.slope)
print('切片 a0: ', lda.intercept)

【実行結果】

線形判別関数の数式表現は次のようになります。

判別得点計算の「関数版」も作りました。ヘルパー関数で定義しています。
引数は ① 判別したいデータ、② 学習データ、③ 変換行列 です。

# 判別得点の算出：関数版
discriminant_score_2d(X_new=X, X_train=X, trans_matrix=lda.e_vecs)

【実行結果】

◆ 可視化
判別得点を可視化しましょう。

## 判別得点の可視化
sns.scatterplot(x=lda_score, y=y,
                hue=np.where(y == 1, '前立腺がん', '前立腺肥大症'),
                palette=['tomato', 'tab:blue'], s=80, alpha=0.9)
plt.axvline(0, color='gray', ls='--', lw=1)
plt.xlabel('判別得点', fontsize=12)
plt.ylabel('グループ', fontsize=12)
plt.ylim(-0.5, 1.5)
plt.yticks([0, 1], labels=[2, 1])
plt.legend(title='実測値', loc='upper left');

【実行結果】

０の垂直線近くにはグループ１の３番目の人が乗っています。
グループ２の 14 番目の人の誤判別が分かります。

3️⃣ 次元削減
線形判別分析は機械学習の文脈で「次元削減」に利用されています。
「主成分分析」の仲間に見えますが、線形判別分析は教師あり学習である点が主成分分析と異なります。

テキストは次元削減について触れていませんが、自作クラスには含めてみました。

モデル構築の際に「n_discriminants=1」と指定したので、１次元に次元削減した結果を算出・表示します。

## データを１次元に次元削減
X_transformed = lda.transform(X).flatten()
X_transformed

【実行結果】

１次元チャートで可視化しましょう。

## 次元削減データの可視化
sns.scatterplot(x=X_transformed, y=y,
                hue=np.where(y == 1, '前立腺がん', '前立腺肥大症'),
                palette=['tomato', 'tab:blue'], s=80, alpha=0.9)
plt.xlabel('変換後データ', fontsize=12)
plt.ylabel('グループ', fontsize=12)
plt.ylim(-0.5, 1.5)
plt.yticks([0, 1], labels=[2, 1])
plt.legend(title='実測値', loc='upper left');

【実行結果】
判別得点とよく似た傾向になっています。

◆ 次元削減の仕組み
次元削減のロジックを確認しましょう。
$${n}$$ 次元に次元削減する場合は、$${n}$$ 列の変換行列を用いて：

とします。
自作クラスでは、次の式で求めた行列の固有ベクトル（固有値の大きい順）を変換行列に使用しています。

$${S_W}$$ はグループ内変動、$${S_B}$$ はグループ間変動であり、テキストでは Section 5.3 以降で登場する統計量です。

「e_vec」属性で変換行列を取り出して、変換計算を試してみましょう。

# １次元の変換ロジック
X @ lda.e_vecs[:, 0]

【実行結果】

こんな感じでシンプルですが…
後で実践する mlxtend ライブラリの次元削減結果と同じ値を得られました！

4️⃣ 予測
「 S さんのデータ」から判別得点を算出して、グループ１／グループ２の判別・分類を予測しましょう。
変数 X_new を使います。

## 16番目の被験者Sさんの予測：クラス１（判別得点がプラス⇒クラス１）
lda.discriminant_score_2d(X_new)

【実行結果】
判別得点がプラスなので、グループ１（前立腺がん）と予測されます。

テキスト p.202 の S さんの判別得点は $${0.111}$$ です。
途中計算の小数点丸め桁数の違いから、差が生じています。

5️⃣ 境界線の可視化
おまたせしました！（してない！？）
２次元チャートに線形判別関数で境界線を引きます！

### 境界線の描画

## 境界線の算出
# 境界線：x軸の値の設定
x_vals = np.linspace(0, 5, 101)
# 傾きと切片の算出
slope, intercept = coef_2d(X, lda.e_vecs)
# 境界線：y軸の値の算出
y_vals = intercept + slope * x_vals

## 描画
# 描画領域の設定
plt.figure(figsize=(6, 6))
# 散布図の描画
sns.scatterplot(data=data1, x='マーカーA', y='マーカーB', s=100,
                hue=np.where(y == 1, '実測：前立腺がん', '実測：前立腺肥大症'),
                palette=['tomato', 'tab:blue'],
                style=np.where(
                    lda_score > 0, '予測：前立腺がん', '予測：前立腺肥大症'))
# 16番目の被験者の散布図の描画
plt.scatter(*X_new.flatten(), marker='*', s=300, color='orange');
# 境界線の描画
plt.plot(x_vals, y_vals, color='tab:red', lw=1, ls='--', zorder=0)
# テキスト
plt.text(x=3.5, y=2, s='グループ$G_1$', fontsize=14)
plt.text(x=1, y=1.2, s='グループ$G_2$', fontsize=14)
# 修飾
plt.legend(title='凡例', bbox_to_anchor=(1, 1))
plt.xlim(0, 5)
plt.ylim(1, 5)
plt.yticks(range(1, 6))
plt.gca().set_aspect('equal')
plt.show()

【実行結果】
青い丸点と赤いｘ点が誤判別です。

次式で示される境界線が赤い点線です。

境界線より上（右）がグループ１（前立腺がん）、下（左）がグループ２（前立腺肥大症）、と判別されます。
S さんの星印は境界線の上（右）なので、グループ１（前立腺がん）と予測されます。

6️⃣ その他
lda インスタンスから取り出せる情報を見ておきます。

◆ グループ内変動 $${S_W}$$

# グループ内変動
lda.SW

【実行結果】

◆ グループ間変動 $${S_B}$$
テキストの線形判別関数はこの「グループ間変動 $${S_W}$$ の最大化」によって導出されています。

# グループ間変動
lda.SB

【実行結果】

◆ $${S_W^{-1} S_B}$$ の固有値

# 固有値
lda.e_vals

【実行結果】

◆ $${S_W^{-1} S_B}$$ の固有ベクトル

# 固有ベクトル（射影行列）
lda.e_vecs

【実行結果】

ひととおり線形判別分析クラスを堪能いたしました。

そういえば、線形判別分析クラスを作ったきっかけをお知らせしていないですね…
実は、テキストの線形判別分析の流派に合致する Python ライブラリを見つけられなかったからです…

（ゆる募）
Pythonの線形判別分析ライブラリをご存じの方、ぜひ教えて下さい！

このあたりで休憩いたしましょう。

mlxtend の線形判別分析

■ mlxtend ライブラリのご紹介
mlxtend は機械学習ライブラリの仲間です。

なのです。

書籍「Python機械学習プログラミング［PyTorch＆scikit-learn編］」の著者の一人、Sebastian Raschka 氏により公開されています。

【 mlxtend 公式サイト】

【書籍リンク】
Python機械学習プログラミング［PyTorch＆scikit-learn編］

私はこの書籍で機械学習をガッツリ学びました。
学びの様子は 🔽 こちらのブログ 🔽 で！（宣伝）

■ mlxtend の線形判別分析 LinearDiscriminantAnalysis の特徴
「次元削減」に特化した作りになっています。
線形判別関数や、判別（分類）予測の機能は無さそうです。

では実装していきましょう！

1️⃣ 線形判別分析の実行
某ライブラリとクラス名「LinearDiscriminantAnalysis」が重なったので、mlxtendの方を「lda_mlxtend」でインポートしました。
線形判別分析のインスタンス生成 ⇒ 線形判別分析の実行の流れです。

## 線形判別分析の実行

# 線形判別分析のインスタンス生成
lda_ml = lda_mlxtend(n_discriminants=1)
lda_ml.fit(X, y)

【実行結果】なし
引数 n_discriminants=1 は、変換（次元削減）時の出力次元数です。
今回は１次元です。

2️⃣ 次元削減
モデル構築の際に「n_discriminants=1」と指定したので、１次元に次元削減した結果を算出・表示します。

## データを１次元に次元削減
X_lda_ml = lda_ml.transform(X)
X_lda_ml.flatten()

【実行結果】
自作クラスと同じ結果になりました。

（参考：自作クラスの出力結果）

変換行列に固有ベクトルを用いて次元削減しているようです。
変換行列（射影行列）を「 .w_ 」で取り出して確認してみましょう。

# LDAの変換は固有ベクトルで行っている
(X @ lda_ml.w_).flatten()

【実行結果】

１次元チャートで可視化しましょう。

## 次元削減データの可視化
sns.scatterplot(x=X_lda_ml[:, 0], y=y,
                hue=np.where(y == 1, '前立腺がん', '前立腺肥大症'),
                palette=['tomato', 'tab:blue'], s=80, alpha=0.9)
plt.xlabel('変換後データ', fontsize=12)
plt.ylabel('グループ', fontsize=12)
plt.ylim(-0.5, 1.5)
plt.yticks([0, 1], labels=[2, 1])
plt.legend(title='実測値', loc='upper left');

【実行結果】

3️⃣ その他
lda_ml インスタンスから取り出せる情報を見ておきます。

◆ 固有値

# 固有値
lda_ml.e_vals_

【実行結果】
自作クラスと違う値になっているのが気になります…
（おそらく内部の計算ロジックが相違しているということでしょう）

◆ 固有ベクトル

# 固有ベクトル
lda_ml.e_vecs_

【実行結果】
自作クラスと同じ値になりました。

（参考：自作クラスの出力結果）

◆ 変換行列（射影行列）

# 射影行列
lda_ml.w_

【実行結果】
固有ベクトルのうち１列目（次元削減後の１次元目）に相当します。

scikit-learn の線形判別分析

■ 特徴
Python の機械学習ライブラリの王道です！
線形判別分析クラス LinearDiscriminantAnalysis() には、次元削減に加えて、判別（分類）予測の機能も備わっています。

🙆‍♂️ いいところ
判別ロジックに「確率」を用いているので、判別の予測時の際、どのグループかの予測に加えて、各グループの確率値も得られます。
高い確率で判別しているのか（判別に自信あり！？）、他のグループをほとんど差がない確率で判別しているのか（判別の自信度は…？）を知ることができます。

⚠️ 注意点
判別ロジックがテキストと異なっています。
そのため、判別結果もテキストと異なります。
scikit-learn は、S さんをグループ２（前立腺肥大症）に分類予測します。

判別ロジックは記事末の付録に書きます。
（書くといっても、ChatGPTの回答を貼り付けただけです…😅）

1️⃣ 線形判別分析の実行
インスタンス生成とモデルの学習の２ステップで線形判別分析を実行します。

## 線形判別分析の実行

# 線形判別分析のインスタンス生成
lda_sk = LinearDiscriminantAnalysis()
# モデルの学習と次元削減変換の実行
X_transformed_sk = lda_sk.fit_transform(X, y)

【実行結果】なし
引数を省略していますが、変換（次元削減）時の出力次元数 n_components はデフォルト計算値の $${1}$$ が設定されています。

2️⃣ 決定関数の算出
決定関数は、どちらのグループに所属するかの予測に用いられます。
scikit-learn の線形判別分析の場合「対数事後確率のグループ間差」が決定関数になります（詳細は記事末の付録で！）。
決定関数$${>0}$$ の場合はグループ１（前立腺がん）、$${<0}$$ の場合はグループ２（前立腺肥大症）を示します。

## 決定関数の算出：モデルの対数事後分布
decision_function_sk = lda_sk.decision_function(X)
decision_function_sk

【実行結果】
被験者No. １～15 の決定関数です。
プラス・マイナスの向きは、判別得点と「ほぼ」同じになっています。
実は被験者No. 3 の人は「グループ１（前立腺肥大症）」と誤判別しています。

（参考：自作クラスの判別得点の出力結果）

◆ 決定関数の確認
lda_sk インスタンスから「傾き」と「切片」の情報を取り出して、決定関数を計算してみましょう。

# 決定関数: sklearnの切片、重みは決定関数に用いられるもの
lda_sk.intercept_ + X @ lda_sk.coef_[0]

【実行結果】

◆ 可視化
決定関数を可視化しましょう。

## 決定関数の可視化
sns.scatterplot(x=decision_function_sk, y=y,
                hue=np.where(y == 1, '前立腺がん', '前立腺肥大症'),
                palette=['tomato', 'tab:blue'], s=80, alpha=0.9)
plt.axvline(0, color='gray', ls='--', lw=1)
plt.xlabel('決定関数', fontsize=12)
plt.ylabel('グループ', fontsize=12)
plt.ylim(-0.5, 1.5)
plt.yticks([0, 1], labels=[2, 1])
plt.legend(title='実測値', loc='upper left');

【実行結果】

グループ１の３番目の人とグループ２の 14 番目の人について、誤判別になっていることが分かります。

◆ 確率の確認
predict_proba() メソッドで「２つのグループの所属確率」を確認しましょう。

# クラス所属確率の確認 lda.predict_proba(X)

# 列名の設定
columns = ['グループ２の確率', 'グループ１の確率']

# クラス所属確率を算出してデータフレーム化
df = pd.DataFrame(lda_sk.predict_proba(X), index=data1['被験者No.'], 
                  columns=columns)

# 確率の大きい方をハイライト表示
(df.style
 .highlight_max(color='bisque', axis=1)
 .format({'グループ２の確率': '{:.3f}', 'グループ１の確率': '{:.3f}'}))

【実行結果】
誤判別の No. ３と 14 は他のデータと比べて、予測したグループの確率値が小さくなっていることが分かります。

3️⃣ 次元削減
１次元に次元削減した結果を算出・表示します。

## データを１次元に次元削減
X_transformed_sk.flatten()

【実行結果】
自作クラスや mlxtend と異なる値になっています。

次元削減データを１次元チャートで可視化しましょう。

## 次元削減データの可視化
sns.scatterplot(x=X_transformed_sk.flatten(), y=y,
                hue=np.where(y == 1, '前立腺がん', '前立腺肥大症'),
                palette=['tomato', 'tab:blue'], s=80, alpha=0.9)
plt.axvline(0, color='gray', ls='--', lw=1)
plt.xlabel('変換行列', fontsize=12)
plt.ylabel('グループ', fontsize=12)
plt.ylim(-0.5, 1.5)
plt.yticks([0, 1], labels=[2, 1])
plt.legend(title='実測値', loc='upper left');

【実行結果】
No. ３のデータは、次元削減後はプラスの値になっているようです。

◆ 次元削減の仕組み
次元削減のロジックを確認しましょう。
中心化したデータに対して変換をかけているようです。
lda_sk インスタンスから変換行列を取り出して、試してみましょう。

この変換行列を用います。

# 射影行列：transformに利用
lda_sk.scalings_

【実行結果】

中心化したデータに変換行列を掛けて変換します。

# transformと一致 ※中心化した変数 @ 変換行列
((X - X.mean(axis=0)) @ lda_sk.scalings_).flatten()

【実行結果】
最初の次元削減結果と同じになりました。

4️⃣ 予測
lda_sk インスタンスに対して「predict」メソッドを適用することで、判別（分類）の予測ができます。
「 S さんのデータ」を予測します。

## 16番目の被験者Sさんの予測 ★テキストの予測結果と異なる
pred = lda_sk.predict(X_new)[0]
print(f"Sさんの予測: {pred} {'前立腺ガン' if pred==1 else '前立腺肥大症'}")

【実行結果】
グループ２（前立腺肥大症）と予測されました。
テキストと結果が異なっています。

クラス所属確率を確認しましょう。

# クラス所属確率 [グループ1, グループ2]
pd.DataFrame(lda_sk.predict_proba(X_new), columns=columns, index=['S']).round(3)

【実行結果】
両方の確率値は近くなっています（判別の自信度は高くないかも？）。

ちなみに S さんデータで決定関数を算出すると…

# 決定関数
lda_sk.decision_function(X_new)

【実行結果】
マイナス値なのでグループ２を予測します。

5️⃣ 境界線の可視化
２次元チャートに境界線を引きます！
決定関数に用いる傾きと切片を利用して、線形判別関数的な関数を暫定作成しています。

### 境界線の描画

## 境界線の算出
# 全体の平均値: 境界線は全体の平均値を通る
mean = X.mean(axis=0)
# 境界線zの重みの取得
w = lda_sk.coef_[0]
# 傾き　wと直交する
b1 = -w[0] / w[1]
# 切片
# b0 = -b1 * mean[0] + mean[1]  # x1,x2の平均値が原点を通る場合
b0 = -lda_sk.intercept_[0] / w[1]  # x1,x2の平均値が原点を通らない場合
# 境界線：x軸の値
x_vals = np.linspace(0, 5, 101)
# 境界線：y軸の値
y_vals = b0 + b1 * x_vals

## 描画
# 描画領域の設定
plt.figure(figsize=(6, 6))
# 散布図の描画
sns.scatterplot(data=data1, x='マーカーA', y='マーカーB', s=100,
                hue=np.where(y == 1, '実測：前立腺がん', '実測：前立腺肥大症'),
                palette=['tomato', 'tab:blue'],
                style=np.where(
                    lda_sk.predict(X)==1, '予測：前立腺がん', '予測：前立腺肥大症'))
# 16番目の被験者の散布図の描画
plt.scatter(*X_new.flatten(), marker='*', s=300, color='orange');
# 境界線の描画
plt.plot(x_vals, y_vals, color='tab:red', lw=1, ls='--', zorder=0)
# テキスト
plt.text(x=3.5, y=2, s='グループ$G_1$', fontsize=14)
plt.text(x=1, y=1.2, s='グループ$G_2$', fontsize=14)
# 修飾
plt.legend(title='凡例', bbox_to_anchor=(1, 1))
plt.xlim(0, 5)
plt.ylim(1, 5)
plt.yticks(range(1, 6))
plt.gca().set_aspect('equal')
plt.show()

【実行結果】
青い丸点（No. 14）と赤いｘ点（No.３）が誤判別です。
S さんはぎりぎりグループ２に判別されます。

（参考：自作クラスの出力結果）

6️⃣ 可視化２：scikit-learn の可視化関数
scikit-learn には 境界線＝決定境界を描画できる DecisionBoundaryDisplay 関数が用意されています。
scikit-learn で２変数の判別（分類）モデルを実装する場合、こちらの関数を使うと境界線が簡単に引けます！

### DecisionBoundaryDisplay利用

## 描画
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 6))
# 境界の領域の色の設定
cmap = colors.ListedColormap(['lightblue', 'peachpuff'])

# 境界の領域の塗りつぶし描画
DecisionBoundaryDisplay.from_estimator(
    lda_sk,
    X,
    response_method='predict_proba',
    plot_method='pcolormesh',
    ax=ax,
    cmap=cmap,
    alpha=0.3,
)

# 境界線の描画
DecisionBoundaryDisplay.from_estimator(
    lda_sk,
    X,
    response_method='predict_proba',
    plot_method='contour',
    ax=ax,
    alpha=1.0,
    cmap='hsv',
    levels=[0.5],
)

# 実測値の散布図の描画
sns.scatterplot(data=data1, x='マーカーA', y='マーカーB', hue='結果',
                palette=['tab:blue', 'tomato'],  s=100)
# Sさんの位置の描画
plt.scatter(x=X_new[0, 0], y=X_new[0, 1], marker='*', s=300, color='tab:orange',
            zorder=10)
# 修飾
plt.gca().set_aspect('equal')
plt.legend(title='実測値');

【実行結果】
どちらのグループに判別されるか「色分け」で区別できるので、直感的に分かりやすいですね！

scikit-learn の線形判別分析は、テキストと異なる結果になるものの、分類予測の確率値を確認できたり、次元削減ができるので、とてもいい感じだと思います。

付録：scikit-learn の線形判別分析を深堀り

「scikit-learn の LinearDiscriminantAnalysis の決定関数」に関するChatGPTの回答をそのまま貼り付けします。

【問１】LDAの決定関数は何か？

【ChatGPTの回答】

「log-posterior の差」というのは、ざっくり言うと

を数値化したものです。以下、もう少しかみ砕いて説明しますね😊

◆

1️⃣ そもそも “Posterior（事後確率）” とは？

(1) 事前確率（Prior）
たとえば、お店に来たお客さんの７割が男性、３割が女性、という情報があったとします。

• 男性である確率：0.7

• 女性である確率：0.3

これが「事前確率（Prior）」。
データを見ずに最初に信じている割合です。

(2) 尤度（Likelihood）
もし「この人は背が高い」「スーツっぽい格好をしている」といった特徴を観察したとき、

• 背が高い人が男性である確率

• スーツを着る人が男性である確率

など、「観察した特徴がそのクラスで出やすいか」を表すものです。

(3) 事後確率（Posterior）
事前確率 × 尤度　をもとにベイズの定理で計算すると、

• たとえば計算結果として「男性である確率は 0.85（85％）」

• 女性である確率は同様に計算して 0.15（15％）

というふうに得られます。

◆

2️⃣ なぜ “Log（対数）” を取るのか？

確率そのもの（たとえば 0.85 や 0.15）を使って比べても結論は同じですが、計算の都合や数値の扱いやすさから「対数」を取ってから比べる場合があります。
理由は主に：

1. 掛け算を足し算に変換したい

   • 事後確率の計算で「いくつもの確率を掛け合わせる」ことが発生します。

   • 連続して掛け算をするよりも、「ログを取って足し算にして計算したほうが安定する」ため。

2. 非常に小さい数字を扱いやすくする

   • 確率をそのまま掛け合わせるととても小さい数字になり、丸め誤差が出やすい。

   • 対数を取ることで、たとえば 0.000001 のような数字が「-13.8…」のような数値になり、
     PC が扱いやすい（ゼロに丸められにくい）メリットがあります。

◆

3️⃣ 「log-posterior の差」とは？

複数クラスがあるとき、LDA（線形判別分析）はクラスごとに

1. 事前確率（Prior）

2. そのクラスでの特徴分布（Likelyhood）

を計算し、最終的に「このデータがクラス $${k}$$ である事後確率（Posterior）を求める」わけですが、scikit-learn の decision_function は 「事後確率を直接返す」のではなく、その「対数を取ったもの（log-posterior）」を返しています。

具体的に言うと、２クラス問題（クラスA と クラスB）では、LDA は

• クラスA の log-posterior： $${\log ⁡P (\text{A} \mid x)}$$

• クラスB の log-posterior： $${\log ⁡P (\text{B} \mid x)}$$

を計算し、その差

を decision_function の値として返します。
（多クラスの場合は、ある「基準クラス」との log-posterior の差を返す形式になっています）

◆

4️⃣ どう解釈するか？

• 差がプラス ⇒ $${\log P(\text{A} \mid x) > \log P(\text{B} \mid x) \Rightarrow P(\text{A} \mid x) > P(\text{B} \mid x)}$$
  → 「クラスA のほうが確率が高い」 → A に分類

• 差がマイナス ⇒ $${P(\text{B} \mid x)}$$ のほうが高い → B に分類

実際に最終分類は np.argmax を使って「最も log-posterior が大きいクラス」を選んでいますが、decision_function を見ると 「A と B のどちらがどれくらい優勢か」 が絶対値でわかります。
（正負だけでなく数値の大きさが「どのくらい差があるか」も示す）

◆

📕 まとめ

1. posterior（事後確率） ＝ 「特徴を見たあとで〈このデータはクラスXである〉確率」

2. その「対数」を取ったものが log-posterior

3. decision_function は、

   • 各クラスの log-posterior を計算し、

   • 基準クラスとの差を返す

4. 正の値なら「このクラスのほうが確率が高い」、負なら「基準クラスのほうが高い」

5. 数字が大きいほど「差が大きい（確信が強い）」と考えられる

これで、「log-posterior の差」が何を表し、なぜ LDA がそれを返すのかがイメージしやすくなるはずです😊

【問２】事後確率の計算に用いる「尤度」は何か？

【ChatGPTの回答】

1️⃣ 尤度（Likelihood）の計算方法

「尤度（likelihood）」というのは、「もしこのデータ点 $${\mathbf{x}}$$ がクラス $${k}$$（例えばクラスA）に属していると仮定したとき、その特徴 $${\mathbf{x}}$$ が得られる確率」をモデル化したものです。
LDA（線形判別分析）では、各クラスのデータが以下のような「多変量正規分布（ガウス分布）」に従うと仮定します。

1. 各クラスの平均ベクトル $${\boldsymbol\mu_k}$$

2. 全クラス共通の共分散行列 $${\Sigma}$$

この仮定のもとで、「クラス $${k}$$ に属するときの尤度」は、次の式のような多変量正規分布の確率密度関数（PDF）として計算されます。

ここで

1. $${(2\pi)^{d/2}|\Sigma|^{1/2}}$$

   • $${(2\pi)^{d/2}}$$ は次元数 $${d}$$ が高くなるにつれて「正規化定数」が大きくなる部分

   • $${|\Sigma|}$$（シグマの行列式）は、共分散行列 $${\Sigma}$$ の「体積」を表す数値。

     • $${\Sigma}$$ が大きなばらつきをもつときは $${|\Sigma|}$$ も大きく、全体的な確率密度が小さくなる。

     • $${\Sigma}$$ がばらつきの小さい方向で小さい固有値をもつときは $${|\Sigma|}$$ が小さくなり、「ピークが高い分布」になる。

2. $${\exp\bigl(-\tfrac{1}{2} (\mathbf{x}-\boldsymbol\mu_k)^\top \Sigma^{-1} (\mathbf{x}-\boldsymbol\mu_k)\bigr)}$$

   • $${\mathbf{x}-\boldsymbol\mu_k}$$ は「データ点がクラス $${k}$$ の平均からどれだけ離れているか」を示す差分ベクトル。

   • $${\Sigma^{-1}}$$（共分散行列の逆行列）をかけることで、各次元のばらつきの違いを考慮した以下の「マハラノビス距離」が現れます。

• つまり「ばらつきの大きい方向では、少しくらい平均から離れても影響が小さく、ばらつきの小さい方向から少し外れると急に確率が下がる」という性質を表します。

• そのマハラノビス距離 $${D_k(\mathbf{x})}$$ が大きいほど（＝平均から外れているほど）、指数部が小さくなり、尤度（確率密度）も小さくなります。

◆

2️⃣ LDA における事後確率への使い方

1. 事前確率 $${\pi_k}$$（クラス $${k}$$ に事前に属してくる確率、例えばデータセット上でクラス $${k}$$ のサンプル割合）

2. 尤度 $${p(\mathbf{x}\mid\text{class }k)}$$（上の式）

をかけ合わせて、ベイズの定理で事後確率を求めます。

• 分母は「すべてのクラスでこの $${\mathbf{x}}$$ が出てくる確率」を合計した正規化項。

• 最終的に $${\displaystyle P(\text{class }k \mid \mathbf{x})}$$ が最も大きいクラスを、「この $${\mathbf{x}}$$ はクラス $${k}$$ に分類される」と判断します。

◆

3️⃣ イメージのまとめ

• 尤度 $${p(\mathbf{x}\mid k)}$$
  → 「この特徴 ($${\mathbf{x}}$$) はクラス (k) の中でどれくらいありそうか？」
  → 多変量正規分布の式で計算する

• 事後確率 $${P(k \mid \mathbf{x})}$$
  → 「$${\mathbf{x}}$$ を見たあとで、このデータはクラス $${k}$$ であるとどれだけ確信できるか？」
  → 尤度×事前確率 をベイズの定理で正規化して求める

• 最後に一番大きい「事後確率」を出すクラスに $${\mathbf{x}}$$ を割り振る

こうして、「特徴 $${\mathbf{x}}$$ がどのクラスから来たデータか」を確率的に判定しているのが LDA の根幹です。
尤度の式の中身を理解すると、「なぜマハラノビス距離で距離を測るのか」「なぜ対数を使うのか」もよりイメージしやすくなります😊

難解です…
テキストの線形判別関数（判別得点）と scikit-learn の決定関数が別物だろう、ということは分かった気がします…

記事の最後はChatGPTが締めくくります。
今回は日々現れる「分かれ道」を振り返りながら。

📘 ChatGPTのひとこと：

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。