「入門はじめての多変量解析」をPythonで写経 Vol.25 ～ ５章「はじめての判別分析」④正答率

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

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

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

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

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

この記事は、正答率 に取り組みます。
判別・分類タスクにおけるモデルの予測精度をはかる指標の一種です。

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

はじめに

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

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

５章 はじめての判別分析

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

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

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

### インポート

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

# 混同行列
from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay

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

正答率・誤判別率

テキストで学んだ判別分析では、「線形判別関数」・「マハラノビス距離」で判別・分類予測をする「モデル」を構築します。
データの個数が限られた標本から 100% 的中するモデルを実現することは難しそうですが（むしろ 100% 精度は怪しい？）、一定レベルの予測精度が出せることには興味があります。

正答率・誤判別率は予測精度をはかる指標の一種です。

データの全個数のうち、正しく判別できた個数の割合が正答率です。
一方で誤判別率は判別を誤った個数の割合です。
つまり、正答率＋誤判別率＝１です。

テキストにならって、例題データ＝学習データを対象にして「線形判別関数」と「マハラノビス距離による判別」で判別予測を行って、正答率を計算しましょう。
※誤判別率は１－正答率で求まりますので省略します！
※バリデーションという概念は取り入れていません。

正答率の算出

1️⃣ データの準備
テキスト p.179 表 5.1.1 のデータを引用いたします。

### マーカー測定結果 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

【実行結果】
結果＝０がグループ２（前立腺肥大症）、結果＝１がグループ１（前立腺がん）です。
標本サイズはグループ１は $${7}$$、グループ２は $${8}$$ です。

データを説明変数と目的変数に分割し、説明変数をグループ別に分けます。

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

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

# グループG1データx_g1とグループG2データx_g2に分割
x_g1 = X[y==1]
x_g2 = X[y==0]

【実行結果】なし

モデルの構築等で利用する統計量を算出します。
全体の変数別平均を算出します。

## 全体の平均値の算出
X_mean = X.mean(axis=0)
print('全体の平均:', X_mean)

【実行結果】

グループ別の平均と分散共分散行列を算出します。

## グループG1の平均・分散共分散行列の算出
# 平均
x_bar_g1 = x_g1.mean(axis=0)
# 分散共分散行列
cov_g1 = np.cov(x_g1, ddof=1, rowvar=False)
# 結果の表示
print('グループ１:')
print(' 平均:')
print(x_bar_g1)
print(' 分散共分散行列:')
print(cov_g1)

## グループG2の平均・分散共分散行列の算出
# 平均値
x_bar_g2 = x_g2.mean(axis=0)
# 分散共分散行列
cov_g2 = np.cov(x_g2, ddof=1, rowvar=False)
# 結果の表示
print('\nグループ２:')
print(' 平均:')
print(x_bar_g2)
print(' 分散共分散行列:')
print(cov_g2)

【実行結果】

線形判別分析とマハラノビス距離に関連する関数を定義します。

線形判別関数の係数算出関数と判別得点算出関数です。

## 線形判別関数の求め方の公式 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

マハラノビス距離（の二乗）算出関数です。

## マハラノビスの距離の二乗を算出する関数  p.211
def mahalanobis_distance_sq(X, x_bar, cov):
    D2 = [(x - x_bar) @ np.linalg.inv(cov) @ (x - x_bar) for x in X]
    return np.array(D2)

2️⃣ データと境界線の可視化
テキスト p.214 図 5.10.1 に相当するチャートで、線形判別関数による境界線とマハラノビス距離による境界線のイメージを膨らませましょう。

### 境界線の描画 p.214 図 5.10.1

## 判別分析による境界線の算出
# 境界線zの係数の取得 ※自作関数ldfを利用
res2 = ldf(x_g1, x_g2)
a0, a1, a2 = res2[0], res2[1], res2[2]
# 傾き: wと直交する
b1 = -a1 / a2
# 切片
b0 = -a0 / a2
# 境界線: x軸の値
x_vals = np.linspace(X[:, 0].min(), X[:, 0].max(), 2)
# 境界線: y軸の値
y_vals = b0 + b1 * x_vals

## マハラノビスの距離による境界線の算出
# X, Yの算出
n = 100
x = np.linspace(1, 5, n)
XX, YY = np.meshgrid(x, x)
XY = np.column_stack([XX.flatten(), YY.flatten()])
# Zの算出 グループG1のマハラノビスの距離とグループG2のマハラノビスの距離の差
ZZ = (mahalanobis_distance_sq(XY, x_bar_g1, cov_g1)
      - mahalanobis_distance_sq(XY, x_bar_g2, cov_g2)).reshape(n, n)

## 描画
# 描画領域の設定
plt.figure(figsize=(6, 6))
# 散布図の描画
sns.scatterplot(data=data1, x='マーカーA', y='マーカーB', s=100,
                hue=np.where(y==1, '実測：グループG1', '実測：グループG2'),
                palette=['tomato', 'tab:blue'])
# 16番目の被験者の散布図の描画
plt.scatter(x=[2.7], y=[3.1], marker='*', s=300, color='orange', label='S')
# 判別分析の境界線の描画
plt.plot(x_vals, y_vals, color='black', lw=1, zorder=0)
plt.text(x=1.8, y=4.5, s='線形判別分析の\n境界線', color='black', fontsize=14)
# マハラノビスの境界線の境界線の描画: 等高線=0=２つのマハラノビスの距離の差を描画する
plt.contour(XX, YY, ZZ, levels=np.arange(0, 1), colors=['tab:red'],
            linestyles=['--'], zorder=0)
plt.text(x=2.3, y=1.1, s='マハラノビス\n距離の境界線', color='tab:red',
         fontsize=14)
# テキスト
plt.text(x=4, y=4.7, s='グループ$G_1$', color='tomato', fontsize=14)
plt.text(x=1.1, y=1.1, s='グループ$G_2$', color='tab:blue', fontsize=14)
# 修飾
plt.ylim(1, 5)
plt.legend(title='凡例', bbox_to_anchor=(1, 1))
plt.show()

【実行結果】
黒実線が「線形判別関数の境界線」、赤破線が「マハラノビス距離の境界線」です。
グループ１・２を区分けする領域が少々異なっています。

3️⃣ 線形判別関数による判別の正答率
テキストにならって、グループ別に正答率を算出します。

◆ グループ１
判別予測の正誤表を作成します。
テキスト p.215 表 5.10.1 に相当します。

### グループG1の判別得点 p.215 表5.10.1
score1 = discriminant_score_2d_with_coefs(x_g1, a0, a1, a2)
g1_df = pd.DataFrame(
    {'判別得点': score1,
     '実測値': y[y==1],
     '予測値': [1 if s>0 else 0 for s in score1]},
     index=data1[data1['結果']==1]['被験者No.'])
g1_df['正誤'] = (g1_df['実測値']==g1_df['予測値']).astype(int)
g1_df.round(3)

【実行結果】
正誤＝１は正答、正誤＝０は誤答です。
計算するまでもなく、正答率は $${1}$$ です。

一応、正答率を計算しましょう。

# 正答率
g1_df['正誤'].mean()

【実行結果】

機械学習ライブラリの scikit-learn はモデル評価指標を算出する便利な関数を多数持っています。
正答率は正解率 accuracy_score で算出します。
引数は、目的変数の正解値、目的変数の予測値、です。

# scikit-learn 利用
accuracy_score(y_true=g1_df['実測値'], y_pred=g1_df['予測値'])

【実行結果】

◆ グループ２
判別予測の正誤表を作成します。
テキスト p.215 表 5.10.2 に相当します。

### グループG2の判別得点 p.215 表5.10.2
score2 = discriminant_score_2d_with_coefs(x_g2, a0, a1, a2)
g2_df = pd.DataFrame(
    {'判別得点': score2,
     '実測値': y[y==0],
     '予測値': [1 if s>0 else 0 for s in score2]},
     index=data1[data1['結果']==0]['被験者No.'])
g2_df['正誤'] = (g2_df['実測値']==g2_df['予測値']).astype(int)
g2_df.round(3)

【実行結果】
誤答が１つあります。

正答率を計算しましょう。

# 正答率
g2_df['正誤'].mean()

【実行結果】
８人中７人が正答ですので、$${7/8＝0.875}$$ です。
誤判別率は $${0.125}$$ です。

scikit-learn の accuracy_score で正答率を算出します。

# scikit-learn 利用
accuracy_score(y_true=g2_df['実測値'], y_pred=g2_df['予測値'])

【実行結果】

◆ データ全体
２つのグループを統合して、データ全体の正答率を算出しましょう。

# 全体の正答率

# ２つのグループのデータフレームを統合
lda_df = pd.concat([g1_df, g2_df], axis=0)
# 正答率の算出
accuracy_score(y_true=lda_df['実測値'], y_pred=lda_df['予測値'])

【実行結果】
全体の正答率は $${0.933}$$ でした。

◆ おまけ：混同行列
グループごとに正解値の個数と予測値の個数をクロス集計したものを「混同行列」と呼ぶようです。
scikit-learn の関数を利用して、混同行列の可視化を試してみましょう。

## 混同行列の可視化

# 混同行列の算出
conf_mtx = confusion_matrix(y_true=lda_df['実測値'], y_pred=lda_df['予測値'])
# 混同行列の可視化の設定
cm = ConfusionMatrixDisplay(conf_mtx, display_labels=['$G_2$', '$G_1$'])
# 混同行列の可視化の実行
cm.plot(cmap='Greens', colorbar=False, text_kw={'fontsize': 30});

【実行結果】
行に「正解値」、列に「予測値」を配したクロス集計表です。

正解値がグループ２の８個のデータのうち、１個をグループ１と予測していることが分かります。

4️⃣ マハラノビス距離による判別の正答率
テキストにならって、グループ別に正答率を算出します。
マハラノビス距離は二乗する定義を使います。

◆ グループ１
判別予測の正誤表を作成します。
判別の指標は「マハラノビス距離の差 $${D_1^2 - D_2^2}$$」です。
指標がマイナス値の場合はグループ１、プラス値の場合はグループ２に判別します。

### マハラノビスの距離による判別の正答率 自作関数利用

## グループG1の距離の差

# グループG1のマハラノビスの距離の差（D₁² - D₂²）の算出
score1_m =   mahalanobis_distance_sq(x_g1, x_bar_g1, cov_g1) \
           - mahalanobis_distance_sq(x_g1, x_bar_g2, cov_g2)

# データフレーム化
g1_m_df = pd.DataFrame(
    {'距離の差': score1_m,
     '実測値': y[y==1],
     '予測値': [1 if s<0 else 0 for s in score1_m]},
     index=data1[data1['結果']==1]['被験者No.'])
g1_m_df['正誤'] = (g1_m_df['実測値']==g1_m_df['予測値']).astype(int)
g1_m_df.round(3)

【実行結果】
すべて正答です。

一応、正答率を計算しましょう。

# 正答率
g1_m_df['正誤'].mean()

【実行結果】

scikit-learn の accuracy_score で正答率を算出します。

# scikit-learn 利用
accuracy_score(y_true=g1_m_df['実測値'], y_pred=g1_m_df['予測値'])

【実行結果】

◆ グループ２
判別予測の正誤表を作成します。

## グループG2の距離の差

# グループG2のマハラノビスの距離の差（D₁² - D₂²）の算出
score2_m =   mahalanobis_distance_sq(x_g2, x_bar_g1, cov_g1) \
           - mahalanobis_distance_sq(x_g2, x_bar_g2, cov_g2)

# データフレーム化
g2_m_df = pd.DataFrame(
    {'距離の差': score2_m,
     '実測値': y[y==0],
     '予測値': [1 if s<0 else 0 for s in score2_m]},
     index=data1[data1['結果']==0]['被験者No.'])
g2_m_df['正誤'] = (g2_m_df['実測値']==g2_m_df['予測値']).astype(int)
g2_m_df.round(3)

【実行結果】
誤答が１つあります。

正答率を計算しましょう。

# 正答率
g2_m_df['正誤'].mean()

【実行結果】
８人中７人が正答ですので、$${7/8＝0.875}$$ です。
誤判別率は $${0.125}$$ です。

scikit-learn の accuracy_score で正答率を算出します。

# scikit-learn 利用
accuracy_score(y_true=g2_m_df['実測値'], y_pred=g2_m_df['予測値'])

【実行結果】

◆ データ全体
２つのグループを統合して、データ全体の正答率を算出しましょう。

# 全体の正答率

# ２つのグループのデータフレームを統合
mah_df = pd.concat([g1_m_df, g2_m_df], axis=0)
# 正答率の算出
accuracy_score(y_true=mah_df['実測値'], y_pred=mah_df['予測値'])

【実行結果】
全体の正答率は $${0.933}$$ でした。

◆ おまけ：混同行列
混同行列を可視化します。

## 混同行列の可視化

# 混同行列の算出
conf_mtx_m = confusion_matrix(y_true=mah_df['実測値'], y_pred=mah_df['予測値'])
# 混同行列の可視化の設定
cm_m = ConfusionMatrixDisplay(conf_mtx_m, display_labels=['$G_2$', '$G_1$'])
# 混同行列の可視化の実行
cm_m.plot(cmap='Greens', colorbar=False, text_kw={'fontsize': 30});

【実行結果】

線形判別関数による判別と同じ結果です。
正解値がグループ２の８個のデータのうち、１個をグループ１と予測していることが分かります。

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

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