こんにちは！データサイエンティストの青木和也（https://twitter.com/kaizen_oni）です！
今回の記事では、統計学の青本「自然科学の統計学」の第9章-演習問題7「線形判別関数による計量分類」についてPythonで解いていきたいと思います。
この問題は、花の特徴量に基づいてどの花である確率が高いかを予測する線形式を使って解いていきます。
問題文
＜線形判別関数による計量分類＞
アイリス・データの各種10ケース計30ケースをランダムに選び出し、それを(9.37)*の線形判別関数で分類してみよ.
分類の過誤（$\prod_i$に属するケースが$\prod_j(i\ne j)$へ分類されること）は生じているか.
*線形判別関数は(9.37)式ではなく(9.42)式と思われるので、(9.42)式を掲載いたします.
$$D_{12}(z)=-3.2456x_1 – 3.3907 x_2 + 7.5530 x_3 + 14.6358 x_4 – 31.5226$$
$$D_{13}(z)=-11.0759x_1 – 19.9161 x_2 + 29.1874 x_3 + 38.4608 x_4 – 18.0933$$
東京大学教養学部統計学教室『自然科学の統計学』(東京大学出版社/2001) 第9章 演習問題 P274
線形判別関数を計算した後は？
$D_{12}(z)$と$D_{13}(z)$の具体的な値を特徴量$x_1, x_2, x_3, x_4$を代入して計算した後は、以下の判別式を使って花の分類を行います。
i ) $D_{12} > 0,~~D_{13}>0$ならば、$\pred_1$（ヴァージニカ）へ
ii ) $D_{12} < 0,~~D_{23} \equiv D_{13} – D_{12}>0$ならば、$\pred_2$（ヴェルシカラー）へ
iii ) $D_{13} < 0,~~D_{23}<0$ならば、$\pred_3$（セトーカ）へそれぞれ分類する
東京大学教養学部統計学教室『自然科学の統計学』(東京大学出版社/2001) 第9章 P272
実装コード
本問題は以下のコードによって解きました。
なお、問題文では各種10サンプル、合計30サンプルとしていますが、今回の検証は全サンプル150件で行っています。
また、アイリス・データはscikit-learnのirisデータを使用しています。
# 基本的なライブラリのインポート
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.datasets import load_iris
# アイリスデータの読み込み
iris = load_iris()
df = pd.DataFrame(iris.data, columns = iris.feature_names)
df['y'] = iris.target
X = df.drop('y', axis = 1)
y = df['y']
def pred_iris(df, target):
# 切片のための定数項
df['const'] = 1
# データフレームを行列に変換
X = df.drop(target, axis = 1)
X = X.to_numpy()
# 線形判別式のベクトルを作成
D12_vec = np.array([-3.2456, -3.3907, 7.5530, 14.6358, -31.5226])
D13_vec = np.array([-11.0759, -19.9161, 29.1874, 38.4608, -18.0933])
# ベクトル->行列
D12_vec = D12_vec.reshape(5,-1)
D13_vec = D13_vec.reshape(5,-1)
# 線形判別式と特徴量から具体的な値を計算
D12 = X @ D12_vec
D13 = X @ D13_vec
df['D12'] = D12
df['D13'] = D13
df['y_pred'] = 0
# 線形判別式の計算結果に応じて、分類する花の種類を決定
# 0: セトーサ
# 1: ヴェルシカラー
# 2: ヴァージニカ
for index, row in df.iterrows():
d12 = row.D12
d13 = row.D13
if d12 > 0 and d13 > 0:
y_pred = 2
elif d12 < 0 and d13 - d12 > 0:
y_pred = 1
else:
y_pred = 0
df.at[index, 'y_pred'] = y_pred
# 観測と予測を比較し、間違っていた場合にフラグ1を立て、合計を計算
df['count'] = df['y'] != df['y_pred']
print('error counts: {}'.format(df['count'].sum()))
display(df)
return df
df = pred_iris(df, 'y')コード解説
今回の問題は以下の3ステップで問題を解いています。
- 特徴量と線形判別関数の係数を行列に変換
- 線形判別関数から具体的な値を計算
- 不等式により、分類すべき花を決定
順を追って解説していきます。
特徴量と線形判別関数の係数を行列に変換
# アイリスデータの読み込み
iris = load_iris()
df = pd.DataFrame(iris.data, columns = iris.feature_names)
df['y'] = iris.target
X = df.drop('y', axis = 1)
y = df['y']
def pred_iris(df, target):
# 切片のための定数項
df['const'] = 1
# データフレームを行列に変換
X = df.drop(target, axis = 1)
X = X.to_numpy()
# 線形判別式のベクトルを作成
D12_vec = np.array([-3.2456, -3.3907, 7.5530, 14.6358, -31.5226])
D13_vec = np.array([-11.0759, -19.9161, 29.1874, 38.4608, -18.0933])
# ベクトル->行列
D12_vec = D12_vec.reshape(5,-1)
D13_vec = D13_vec.reshape(5,-1)
pred_iris関数の中で、アイリス・データの特徴量を抽出し、以下2つの処理を行っています。
- 定数項1を追加する
- データフレームをNumpyの行列に変換する
これらの操作は後々に行列積を計算するための準備です。
また、線形判別関数についても係数をそれぞれベクトル化し、行列積を計算できるように(5,1)の行列になるように調整します。
線形判別関数から具体的な値を計算
# 線形判別式と特徴量から具体的な値を計算
D12 = X @ D12_vec
D13 = X @ D13_vec
df['D12'] = D12
df['D13'] = D13特徴量と線形判別関数の係数の行列積を取ることによって、具体的な値を計算します。
計算した後はベクトルが出力されるので、データフレームに格納します。
不等式により、分類すべき花を決定
df['y_pred'] = 0
# 線形判別式の計算結果に応じて、分類する花の種類を決定
# 0: セトーサ
# 1: ヴェルシカラー
# 2: ヴァージニカ
for index, row in df.iterrows():
d12 = row.D12
d13 = row.D13
if d12 > 0 and d13 > 0:
y_pred = 2
elif d12 < 0 and d13 - d12 > 0:
y_pred = 1
else:
y_pred = 0
df.at[index, 'y_pred'] = y_pred
# 観測と予測を比較し、間違っていた場合にフラグ1を立て、合計を計算
df['count'] = df['y'] != df['y_pred']
print('error counts: {}'.format(df['count'].sum()))
display(df)
return df
df = pred_iris(df, 'y')先ほど計算した$D_{12}(z)$と$D_{13}(z)$の値と以下判別式を使って、花の分類を行う
i ) $D_{12} > 0,~~D_{13}>0$ならば、$\pred_1$（ヴァージニカ）へ
ii ) $D_{12} < 0,~~D_{23} \equiv D_{13} – D_{12}>0$ならば、$\pred_2$（ヴェルシカラー）へ
iii ) $D_{13} < 0,~~D_{23}<0$ならば、$\pred_3$（セトーカ）へそれぞれ分類する
東京大学教養学部統計学教室『自然科学の統計学』(東京大学出版社/2001) 第9章 P272
得られたエラー数を見てみると、3つほど誤分類（分類の過誤）が観測されているものの、$3/150 = 2%$と発生確率はかなり低いことがわかる
まとめ
今回の記事では、統計学の青本「自然科学の統計学」の第9章-演習問題7「線形判別関数による計量分類」についてPythonで解いていきました。
本文中の線形判別関数を使って解いていきましたが、多変量正規分布の仮定を置いたときになぜあのような線形判別関数が導けるのかについてはブログ主も把握できていないので、多変量解析について学習を進めたのちに本記事は更新させていただきたいと思います！
コメント