「入門はじめての多変量解析」をPythonで写経 Vol.23 ～ ５章「はじめての判別分析」②線形判別関数・エピローグ

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

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

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

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

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

この記事ではテキストに沿って線形判別関数を求めます。

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

はじめに

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

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

５章 はじめての判別分析

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

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

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

### インポート

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

# LaTeX表示
from IPython.display import Math

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

線形判別関数の算出方法

■ 算出方法と進め方
この記事は、２変数データに関する線形判別関数 $${z = a_1x_1+a_2x_2+a_0}$$ の係数 $${a_0, a_1, a_2}$$ の算出に取り組みます。

テキストは ２つの係数の算出方法を掲載しています。
1️⃣ 判別得点の全変動の中でグループ間変動を最大化
2️⃣ 線形判別関数の求め方の公式

この記事は、次のように実践しようと思います。
1️⃣ はテキストの量多めの数式展開を Python で書き換えていきます。
2️⃣ はテキストの公式を Python 関数にします。

■ 分析の準備
テキストの例題データをお借りします。

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

data1 = pd.DataFrame(
    {'被験者No.': range(1, 16),
     'マーカーA': [3.4, 3.9, 2.2, 3.5, 4.1, 3.7, 2.8,
                  1.4, 2.4, 2.8, 1.7, 2.3, 1.9, 2.7, 1.3],
     'マーカーB': [2.9, 2.4, 3.8, 4.8, 3.2, 4.1, 4.2,
                  3.5, 2.6, 2.3, 2.6, 1.6, 2.1, 3.5, 1.9],
     '結果': np.hstack([np.ones(7), np.zeros(8)]).astype(int)})
data1

【実行結果】
15人の測定結果です。
結果＝1がグループ１（前立腺がん）、結果＝0がグループ２（前立腺肥大症）です。

学習データと予測データを設定します。
numpy 配列に変換しています。

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

# 学習データの説明変数Xと目的変数y
X = data1[['マーカーA', 'マーカーB']].values
y = data1['結果'].values

# 予測データ（Sさんのデータ） shape(n, d)
X_new = np.array([[2.7, 3.1]])

【実行結果】なし

線形判別分析クラスとヘルパー関数を定義します。
前回記事で紹介したものです！

## 線形判別分析クラスの定義
# 参考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]

    # 変換関数（固有ベクトルを用いる次元削減） ※mlxtedと一致する
    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]
        # 傾きの算出
        self.slope = w[0] / w[1]
        # 切片の算出
        self.intercept = - (self.slope * self.means[0] + self.means[1])
        # 判別得点の算出：線形判別関数 z = a1*x1 + x2 + a0 で計算
        z = self.slope * X_new[:, 0] + X_new[:, 1] + self.intercept
        # 戻り値：判別得点z
        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

判別得点の全変動の中でグループ間変動を最大化

判別得点は、線形判別関数 $${z = a_1x_1+a_2x_2+a_0}$$ の $${x_1, x_2}$$ にデータの具体的な値を当てはめて算出される統計量です。

■ ゴール
テキストが目指す「最も良い判別」の条件

ために、判別得点に関する

を算出します。
「変動」は徐々に明らかになります！

■ 判別得点のイメージ
判別得点的なものを図示しましょう。

線形判別分析を実行して、ひとまず線形判別関数を求めます。
答え先出し方式です。

## 線形判別分析の実行

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

【実行結果】なし

判別得点的なものを可視化します。
テキスト p.185 図 5.3.1 に類するチャートです。

### 境界線の描画 p.185 図5.3.1 ＋ p.202 図5.6.1

## 境界線の算出
# 境界線：x軸の値
x_vals = np.linspace(X[:, 0].min(), X[:, 1].max(), 2)
# 傾きと切片の算出
slope, intercept = coef_2d(X, lda.e_vecs)
# 境界線：y軸の値の算出
y_vals = intercept + slope * x_vals

## 境界線への直交射影点の算出
# データの中心点（平均点）の算出
X_mean = X.mean(axis=0)
# 境界線への射影行列の算出
Base = np.array([[1, slope]]).T
Proj_mtx = Base @ np.linalg.inv(Base.T @ Base) @ Base.T
# データを境界線へ直交射影
X_proj = (X - X_mean) @ Proj_mtx.T + X_mean

## その他の設定
# 予測値の算出
y_pred = (lda.discriminant_score_2d(X) > 0).astype(int)
# データ点の色と形状の設定に用いる値
obs_hue = np.where(y==1, '実測：$G_1$', '実測：$G_2$')
pred_style = np.where(y_pred==1, '予測：$G_1$', '予測：$G_2$')


## 描画
# 描画領域の設定
plt.figure(figsize=(6, 6))

# 散布図の描画
sns.scatterplot(data=data1, x='マーカーA', y='マーカーB', s=100,
                hue=obs_hue, palette=['tomato', 'tab:blue'],
                style=pred_style, zorder=10)

# 16番目の被験者の散布図の描画
plt.scatter(*X_new.flatten(), marker='*', s=100, color='orange', label='S')

# 境界線の描画
plt.plot(x_vals, y_vals, color='tab:red', lw=1, ls='--', zorder=0)
plt.text(x=1, y=4.5, s='境界線', color='tab:red', fontsize=14)

# 直交射影点の描画
plt.plot(*X_proj.T, 'o', color='tab:red', ms=3.5)
for dot1, dot2 in zip(X, X_proj):
    plt.plot(*np.column_stack([dot1, dot2]), lw=1, ls=':', color='black',
             zorder=0)

# データの中心点の描画
plt.plot(*X_mean, 'D', color='black', alpha=0.7, label=f'中心{X_mean.round(3)}',
         zorder=0)

# テキスト
plt.text(x=3.5, y=2, s='グループ$G_1$', color='tomato', fontsize=14)
plt.text(x=1, y=1.2, s='グループ$G_2$', color='tab:blue', 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()

【実行結果】

赤い点線が境界線 $${0=a_1x_1+a_2x_2+a_0}$$ です。
丸やバツ点のデータ点から境界線に向かって垂線を引いています。
この垂線が判別得点的なものになります。

正確には…テキストのヘッセの標準形を用いて、
データ点 $${p, q}$$ と境界線の距離は

となります。
判別得点は $${a_1p + a_2q + a_0}$$ です。
ですので、先ほどのチャートの垂線の長さに $${\boldsymbol{\sqrt{a_1^2 + a_2^2}}}$$ を掛けると判別得点になります。

垂線の長さ×$${\sqrt{a_1^2 + a_2^2}}$$ = 判別得点を検算してみます。
先に答えとしての判別得点を算出します。

# 判別得点
disc_score = lda.discriminant_score_2d(X)
disc_score

【実行結果】

続いて、垂線の長さ×$${\sqrt{a_1^2 + a_2^2}}$$ = 判別得点 を計算します。

# 判別得点の図形的意味 ※ヘッセの標準形 * √(a₁² + a₂²) * 符号
a1, a2 = slope, 1
X_subs = X - X_proj
np.linalg.norm(X_subs, axis=1) * np.sqrt(a1**2 + a2**2) * np.sign(X_subs[:, 0])

【実行結果】
一致しました。検算完了です。

■ 判別得点と変動
判別得点に関する変動は「個々のデータの判別得点と判別得点の平均との差の二乗和」です。
変動は「全変動 $${S_T}$$」「グループ間変動 $${S_B}$$」「グループ内変動 $${S_W}$$」で表現されます。

アルファベットの略はおそらく…

• S ：二乗和を示す Sum of Squares の どちらかの S

• T ：Total（総）

• B ：Between groups（グループ間）

• W ：Within groups（グループ内）

グループ数が２つのケースを対象にして変動を確認します。

① 全変動 $${S_T}$$
判別得点と判別得点の全平均との差の二乗和です。
テキストの数式は、グループごとに二乗和を算出することを意図しているようです。
「判別得点の全平均 $${\bar{z}}$$」、「グループ１の判別得点 $${z_i^{(1)}}$$」（$${i=1, \cdots, N_1}$$）、「グループ２の判別得点 $${z_i^{(2)}}$$」（$${i=1, \cdots, N_2}$$）を用いて、全変動は

です。

② グループ間変動 $${S_B}$$
「各グループの判別得点の平均 $${\bar{z}^{(1)}, \bar{z}^{(2)}}$$」と「判別得点の全平均 $${\bar{z}}$$」の差の二乗を、それぞれグループごとの標本サイズ $${N_1, N_2}$$ 分足し上げるものです。

③ グループ内変動 $${S_W}$$
「各グループの判別得点」と「各グループの判別得点の平均」の差の二乗和です。

④ ３つの変動の関係
ずばり、「全変動＝グループ間変動＋グループ内変動」です。

⑤ 条件の定式化
ゴールに掲げた

について、テキストは変動の比 $${F}$$ を用いて、

で定式化しています。
$${F(a_1, a_2)}$$ が最大となる $${a_1, a_2}$$ を求めに行きます！

線形判別関数の算出 p.186～

テキスト p.186 ～ 199 の計算を Python で一気に実践します！
sympy ライブラリを利用します。

0️⃣ 準備
係数 $${a_0, a_1, a_2}$$ を設定します。

### sympyで計算を進める

# 変数の定義
a0, a1, a2 =  sympy.symbols('a0 a1 a2')

【実行結果】なし

1️⃣ 判別得点の３つの平均の算出
グループ１の判別得点の平均 $${\bar{z}^{(1)}}$$ p.186

# z1_meanの算出
z1_coef = X[y==1].mean(axis=0)
z1_mean = z1_coef[0] * a1 + z1_coef[1] * a2 + a0
display(Math(f'\\bar{{z}}^{{(1)}} = {sympy.latex(z1_mean)}'))

【実行結果】

【TeX 表示の補足説明】
Math で数式を TeX 表示します。
sympy の数式は sympy.latex() で囲って上げて、TeX構文に変換します。

グループ２の判別得点の平均 $${\bar{z}^{(2)}}$$ p.187

# z1_meanの算出
z2_coef = X[y==0].mean(axis=0)
z2_mean = z2_coef[0] * a1 + z2_coef[1] * a2 + a0
display(Math(f'\\bar{{z}}^{{(2)}} = {sympy.latex(z2_mean)}'))

【実行結果】

全体の判別得点の平均 $${\bar{z}}$$ p.187

# z_meanの係数
z_coef = X.mean(axis=0)
z_mean = z_coef[0] * a1 + z_coef[1] * a2 + a0
display(Math(f'\\bar{{z}} = {sympy.latex(z_mean)}'))

【実行結果】

2️⃣ ３つの変動の算出
◆ グループ１のグループ内変動 p.192
リスト内包表記を用いて $${x_1, x_2}$$ の１つ１つの平均値との差の二乗を計算します。

### グループ内変動SWの計算 p.192

# グループG1内の変動の算出
SW_G1 = sum([((a1*x1 + a2*x2 + a0) - z1_mean)**2 for x1, x2 in X[y==1]])
print('【グループG1内の変動】')
display(SW_G1.expand())

【実行結果】

◆ グループ２のグループ内変動 p.193

# グループG2内の変動
SW_G2 = sum([((a1*x1 + a2*x2 + a0) - z2_mean)**2 for x1, x2 in X[y==0]])
print('【グループG2内の変動】')
display(SW_G2.expand())

【実行結果】

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

# ２つのグループ内変動を合計する
SW = SW_G1 + SW_G2
display(Math(f'S_W = {sympy.latex(SW.expand())}'))

【実行結果】

◆ グループ間変動 $${S_B}$$ p.194

### グループ間変動SBの計算 p.194
SB =   ((z1_coef[0]*a1 + z1_coef[1]*a2 + a0) - z_mean)**2 * sum(y==1) \
     + ((z2_coef[0]*a1 + z2_coef[1]*a2 + a0) - z_mean)**2 * sum(y==0)

display(Math(f'S_B = {sympy.latex(SB.expand())}'))

【実行結果】

◆ 全変動 $${S_T}$$ p.194

### 全変動STの計算 p.194
ST = SW + SB
display(Math(f'S_T = {sympy.latex(ST.expand())}'))

【実行結果】

3️⃣ 最適化対象の関数 $${F(a_1, a_2)}$$ p.196

# F(a1, a2)を算出 p.196
F = (SB / ST).simplify()
Math(f'F(a_1, a_2) = {sympy.latex(F)}')

【実行結果】

4️⃣ $${F(a_1, a_2)}$$ を $${a_1, a_2}$$ で偏微分
◆ $${a_1}$$ の偏微分
sympy.diff( 関数, 変数 ) で関数について変数に関する偏微分を計算します。

# Fをa1で偏微分 p.197
expr1 = sympy.diff(F, a1)
expr1.simplify()

【実行結果】

あれあれ？雲行きが怪しくなってきました…

◆ $${a_2}$$ の偏微分

# Fをa2で偏微分 p.197
expr2 = sympy.diff(F, a2)
expr2.simplify()

【実行結果】

・・・・・（大丈夫かな…）

5️⃣ 偏微分の結果を０とおいて連立方程式を解く

# 連立方程式を解く・・・解けない・・・
sympy.solve([expr1, expr2])

【実行結果】

解けませんでした～😭

作戦を変更します！
ChatGPTに解き方を訊きます！

線形判別関数の算出・リベンジ

ここからは ChatGPT の指南にしたがって Python コードを逐次書きます。

1️⃣ 一般には「レイリー商」の極値問題とみなす
３番の $${F(a_1, a_2)}$$ に戻ります。

$$
F = \frac{6.39629761904762 a_1^2 + 10.9077380952381 a_1a_2 + 4.65029761904762 a_2^2}
{11.3293333333333 a_1^2 + 8.76666666666667 a_1a_2 + 12.2133333333333 a_2^2}
$$

は、二次形式

のかたちになっています。

二次形式になっていることを、分子、分母を計算して、確かめましょう。

◆ 行列 $${M}$$ の設定

# 行列 M の設定
M = sympy.Matrix([[SB.expand().coeff(a1, 2), SB.expand().coeff(a1*a2, 1)/2],
                  [SB.expand().coeff(a1*a2, 1)/2, SB.expand().coeff(a2, 2)]])
display(Math(f'M = {sympy.latex(M)}'))

【実行結果】

◆ 行列 $${N}$$ の設定

# 行列 N の設定
N = sympy.Matrix([[ST.expand().coeff(a1, 2), ST.expand().coeff(a1*a2, 1)/2],
                  [ST.expand().coeff(a1*a2, 1)/2, ST.expand().coeff(a2, 2)]])
display(Math(f'N = {sympy.latex(N)}'))

【実行結果】

◆ ベクトル $${a}$$ の設定

# ベクトル a の設定
a = sympy.Matrix([[a1], [a2]])
display(Math(f'a = {sympy.latex(a)}'))

【実行結果】

◆ 分子 $${a^{\top} M a}$$ の算出

# a.T @ M @ a の確認
display(Math(f'a^{{\\top}} M a = {sympy.latex((a.T @ M @ a).expand())}'))

【実行結果】
$${F}$$ の分子と一致しました！

◆ 分母 $${a^{\top} N a}$$ の算出

# a.T @ N @ a の確認
display(Math(f'a^{{\\top}} N a = {sympy.latex((a.T @ N @ a).expand())}'))

【実行結果】
$${F}$$ の分母と一致しました！

確かに $${F}$$ はChatGPTのとおり、二次形式になっています。

2️⃣ 解法～固有値と固有ベクトル
$${\nabla_a F = 0}$$ の条件は、古典的に「一般化固有値問題」

の解（固有ベクトル）を求めることと同等です。
非自明解 $${a \neq [0,0]^\top}$$ のためには

を満たす $${\lambda}$$ を探し、その $${\lambda}$$ に対応するベクトル $${a}$$ を得ます。

ということで…

3️⃣ 固有値 $${\lambda}$$ を求める

を解いて、固有値 $${\lambda}$$ を求めます。
sympy.det() で行列式を計算します。
その後、sympy.solve() で行列式の方程式の解を求めます。

# 数値計算で λ を求める

# 変数の追加定義
λ = sympy.symbols('λ')

# 行列式を定義
det_MN = sympy.det(M - λ*N)
display(Math(f'\det (M - \lambda N) = {sympy.latex(det_MN)}'))

# 行列式を解く
solve_λ = sympy.solve(det_MN)
display(Math(f'\lambda_1, \lambda_2 = {sympy.latex(solve_λ)}'))

# λ₂の取得
λ2 = solve_λ[1]
display(Math(f'\lambda_2 = {sympy.latex(λ2)}'))

【実行結果】

数値計算の結果はおおよそ

となりました。実際のところ $${\lambda_1}$$ はほぼゼロ（行列値の丸め誤差）とみなしますので、興味があるのは $${\lambda_2=0.6965080\cdots}$$ のほうです。

4️⃣ $${\lambda_2}$$ に対応する固有ベクトルを求める

の左辺に $${M, N, \lambda_2}$$ の値を代入すると

# (M - λ₂N)a の設定
expr_λ2 = (M - λ2*N) @ sympy.Matrix([[a1], [a2]])
display(Math(f'(M - \lambda_2 N) a = {sympy.latex(expr_λ2)}'))

【実行結果】

です。例えば上段（１行目）の式は

となります。これを整理すると

$${a_2 = 1}$$ とおいて、係数 $${a_1, a_2}$$ は

となります。

Pythonも追いかけます。
１行目の式に「subs」メソッドで $${a_2 = 1}$$ を代入した上で、sympy.solve() で方程式の解を求めます。

# １行目についてa2=1を代入したうえで、a1を解く
solve_a1 = sympy.solve(expr_λ2.subs(a2, 1)[0])[0]
display(Math(f'a_1 = {sympy.latex(solve_a1)},\quad a_2 = 1'))

【実行結果】

係数のうち $${a_1, a_2}$$ が求まりました！

5️⃣ $${a_0}$$ を求めてゴール！

切片 $${a_0}$$ は、

で算出できます（p.199）

# 線形判別関数のx1の係数slopeと切片interceptの算出
slope = float(solve_a1)
intercept = X_mean @ -np.array([slope, 1])
display(Math(f'a_1 = {slope},\quad a_2 = 1, \quad a_0 = {intercept}'))

【実行結果】

線形判別関数のすべての係数が求まりました！

お疲れ様でした！

線形判別関数の求め方の公式 p.200～

テキストは p.200 ～ 201 で２変数・２グループの場合の「線形判別関数の求め方の公式」をまとめています。

📊 グループごとの分散共分散行列：

を用いて、
📊 プールされた分散共分散行列 $${S}$$：

を算出し、以下の公式で係数 $${a_1, a_2}$$ を求めます。

📊 線形判別関数の係数 係数 $${a_1, a_2}$$

全データの変数ごとの平均 $${\bar{x}_1, \bar{x}_2}$$ を用いて、係数 $${a_0}$$ を求めます。

📊 線形判別関数の係数 $${a_0}$$

これらの公式を関数化しましょう。
算出した係数を用いる判別得点算出関数もあわせて定義します。

### 線形判別関数の求め方の公式 p.200～201

# 線形判別関数の係数を算出する関数（２グループ、２変数）
def ldf(g1, g2):
    
    ## 設定と準備
    # グループ1,グループ2の標本サイズ
    N1, N2 = len(g1), len(g2)
    # グループ1,グループ2の標本平均
    x_bar_g1, x_bar_g2 = g1.mean(axis=0), g2.mean(axis=0)
    # 同じ変数の標本平均の差
    x_bar_diff = x_bar_g1 - x_bar_g2
    # 全データの変数ごとの標本平均
    x_bar = np.vstack([g1, g2]).mean(axis=0)
    
    ## プールされた分散共分散行列の算出
    # グループ1,グループ2の分散共分散の算出
    cov1 = np.cov(g1.T, ddof=1)
    cov2 = np.cov(g2.T, ddof=1)
    # プールされた分散共分散行列の算出
    S = ((N1 - 1) * cov1 + (N2 - 1) * cov2) / (N1 + N2 - 2)
    
    ## 係数a0,a1,a2の算出
    # a1,a2の算出（S=x_bar_diffの連立方程式を解く）
    solve_a1a2 = np.linalg.solve(S, x_bar_diff)
    # a0の算出
    a0 = - x_bar @ solve_a1a2
    
    ## 戻り値: 線形判別関数zの係数a0, a1, a2, プールされた分散共分散行列S
    return a0, solve_a1a2[0], solve_a1a2[1], S

# 判別得点算出関数（上記関数で求めた係数を利用）
def discriminant_score_2d_with_coefs(X, a0, a1, a2):
    z = a1 / a2 * X[:, 0] + X[:, 1] + a0 / a2
    return z

【実行結果】なし

この関数を用いて、例題データの係数を算出します。

### 線形判別関数の係数を算出
a0, a1, a2, S = ldf(X[y==1], X[y==0])
print(a0, a1, a2)

【実行結果】

この計算結果は、テキスト p.200 のうさぎさんの吹き出しコメントに含まれる「係数」と一致しています。

$${a_2 = 1}$$ と置く場合は

# a2=1 とする場合
print(a0/a2, a1/a2, a2/a2)

【実行結果】

プールした分散共分散行列を表示しましょう。
テキスト p.201 に掲載されています。

# プールされた分散共分散行列 p.201
S

【実行結果】

あらためて判別得点を読む

「線形判別関数の求め方の公式」の関数を用いて、あらためて判別得点をまとめます。

グループ１の判別得点です。
テキスト p.202 表 5.6.1 に相当します。

### グループG1の判別得点 p.202 表5.6.1
pd.DataFrame(discriminant_score_2d_with_coefs(X[y==1], a0, a1, a2),
             index=data1[y==1]['被験者No.'], columns=['判別得点']).round(3)

【実行結果】
すべて正の値になっています。

続いてグループ２の判別得点です。
テキスト p.202 表 5.6.2 に相当します。

### グループG2の判別得点 p.202 表5.6.2
df = pd.DataFrame(discriminant_score_2d_with_coefs(X[y==0], a0, a1, a2),
                  index=data1[y==0]['被験者No.'], columns=['判別得点'])
(df.style
.set_properties(
    **{'background-color': 'yellow'}, 
    subset=pd.IndexSlice[df[df['判別得点'] > 0].idxmax(), :])
.format({'判別得点': '{:.3f}'}))

【実行結果】
ほぼ負の値になっていますが、No.14 は正の値（誤判別）となっています。

線形判別関数による境界線を図示します。
テキスト p.202 図 5.6.1 に相当します。

### 境界線の描画 p.202 図5.6.1

## 境界線の算出
# 境界線：x軸の値
x_vals = np.linspace(X[:, 0].min(), X[:, 1].max(), 2)
# 傾きと切片の算出
slope, intercept = coef_2d(X, lda.e_vecs)
# 境界線：y軸の値の算出
y_vals = intercept + slope * x_vals

## その他の設定
# 予測値の算出
y_pred = (discriminant_score_2d_with_coefs(X, a0, a1, a2) > 0).astype(int)
# データ点の色と形状の設定に用いる値
obs_hue = np.where(y==1, '実測：$G_1$', '実測：$G_2$')
pred_style = np.where(y_pred==1, '予測：$G_1$', '予測：$G_2$')


## 描画
# 描画領域の設定
plt.figure(figsize=(6, 6))

# 散布図の描画
sns.scatterplot(data=data1, x='マーカーA', y='マーカーB', s=100,
                hue=obs_hue, palette=['tomato', 'tab:blue'],
                style=pred_style, zorder=10)

# 16番目の被験者の散布図の描画
plt.scatter(*X_new.flatten(), marker='*', s=100, color='orange', label='S')

# 境界線の描画
plt.plot(x_vals, y_vals, color='tab:red', lw=1, ls='--', zorder=0)
plt.text(x=1, y=4.5, s='境界線', color='tab:red', fontsize=14)

# テキスト
plt.text(x=3.5, y=2, s='グループ$G_1$', color='tomato', fontsize=14)
plt.text(x=1, y=1.2, s='グループ$G_2$', color='tab:blue', fontsize=14)
# 修飾
plt.legend(title='凡例', bbox_to_anchor=(1.27, 1))
plt.xlim(0, 5)
plt.ylim(1, 5)
plt.yticks(range(1, 6))
plt.gca().set_aspect('equal')
plt.show()

【実行結果】

S さんの判別得点を算出します。

# Sさんの判別得点
discriminant_score_2d_with_coefs(X_new, a0, a1, a2)

【実行結果】
正の値ですので、グループ１に判別されました。
上のチャートでも境界線の上側＝グループ１にプロットされています。

前回記事・今回記事で線形判別分析に取り組みました。
面白かったですね！

ところで…

📢 判別分析の深堀り情報のご紹介 📢
線形判別分析を数式面で詳細に解説する X の投稿を紹介いたします！
参考になると思いますので、ぜひご一読下さい。

記事の最後は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の教科書です。
よかったらぜひ、お試しくださいませ。

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