判別分析／Python

判別分析

Fisher’s Linear Discriminant Function

各個体から変量のデータが得られていて、これらの個体が個のグループに分けられているとき、グループ間の変動が最大（極大）となる変量の線形結合を求める判別分析のPythonスクリプトを作成した。本ウェブサイトにおける判別分析の説明は、このドキュメントに用意したが、概略は以下の通りである。また、サンプルプログラムとサンプルデータファイルは、まとめて圧縮ファイルとした。サンプルプログラムとその使用例は、このウェブサイトの後半に示す。

第グループの番目の個体のデータを

と表す。は第グループの番目の個体の番目の変量の値を表す。

データ全体の平均を

とおく。ここで、であり、は第グループのデータ（個体）の総数を表す。なお、プログラムではを標準化した値を求めて、この標準化した値に対して判別分析を行っている。

　第グループの平均を以下のようにおく。

データを中心化した（平均を０とした）ものを

とおく。このとき、のデータ全体にわたる平均は

であり、グループにおける平均は

となる。

　判別関数をの変量の線形結合

として構成し、のグループ間の変動とグループ内の変動との比が最大になるようにを決める。グループ間の変動は以下のように算出される。

ここで、

である。

　グループ内変動は以下のように算出される。

ここで、

である。

　級間変動と級内変動の比を

とおく。の最大化を

の条件の下で行う。このために次の目的関数をおく。

　の偏導関数は次式で与えられる。

　上式の値を０とおいて次式を得る。

（１）　

　（１）式を満たすを求める。まず、の固有値分解を考える。

（２）　

　（２）式を（１）式に代入すると以下のようになる。

したがって

（３）　

　ここで、

が成り立っている。すなわち、とは実対称行列の固有値と固有ベクトルになっている。

　（３）式を満たすとは（１）式より次式を満たす。

　したがって、

が成り立つ。すなわち、固有値は群間変動量と群内変動量の比の値になっている。

　0でない固有値の数をとおき、個の固有値は

と並べられているとする。

　番目の固有値に対する固有ベクトルによって構成される判別関数の分散は以下のように計算される。

とおけば、

となる。ここで、

はデータの全変動量を表す。

　したがって、番目の判別関数を平均が０、分散が１になるように標準化したもの（因子、factor）は以下のように与えられる。

　番目までの判別関数を採用したとき、残りの個の判別関数による判別が有意であることの検定は次式により行うことができる。

（４）　、　

ここで、

であり、はにおける変数（変量）の数、はグループの数、はデータの総数である。

　のときの（４）式は、グループ間の差を１元配置の多変量分散分析で検定するものとなっている。

サンプルプログラムと使用例

プログラムは、以下の通りであり、必要なサブプログラムとサンプルデータはまとめて圧縮ファイルとした。実行例は、プログラムリストの後に示す。

"""

Linear Discriminant Analysis

Y. Okamoto, 2019

"""

import numpy as np

import csv

from scipy import stats as ss

import seaborn as sb

import matplotlib.pyplot as plt

fin_name = input('Data file (*.csv) = ')

fout_name = input('Output file (*.txt) = ')

fout = open(fout_name, 'w')

fout.write('Data file = {}\n'.format(fin_name))

with open(fin_name) as fin:

data = [v for v in csv.reader(fin)]

print('data...\n', data)

n_dim = len(data[0]) - 2

VarID = []

for i in range(n_dim):

VarID.append(data[0][i + 2])

print('VarID...\n', VarID)

n = len(data) - 1

CaseID = []

y = []

X = []

for i in range(n):

CaseID.append(data[i + 1][0])

y.append(int(data[i + 1][1]))

X.append([])

for j in range(n_dim):

X[i].append(float(data[i + 1][j + 2]))

y = np.array(y)

X = np.array(X)

fout.write('\n')

for i in range(n):

print(CaseID[i], y[i], X[i])

fout.write('{0} {1} {2}\n'.format(CaseID[i], y[i], X[i]))

cnts = np.bincount(y)

n_groups = len(cnts) - 1

print('n_groups = ', n_groups)

fout.write('\nn_groups = {}\n'.format(n_groups))

n_i = cnts[1:]

print('n_i =\n', n_i)

if n != np.sum(n_i):

print('Counting error')

raise Exception()

Z = ss.zscore(X, axis = 0)

m_g = np.mean(Z, axis = 0)

print('m_g =\n', m_g)

m_i = np.empty((n_groups,n_dim))

for j in range(n_groups):

m_i[j] = np.mean(Z[y == j + 1], axis = 0)

fout.write('\n')

for j in range(n_groups):

print('m_i[{}] = \n{}'.format(j + 1, m_i[j]))

fout.write('\nm_i[{}] = \n{}'.format(j + 1, m_i[j]))

fout.write('\n')

# The between-group scatter matrix

A = np.zeros((n_dim, n_dim))

for i in range(n_groups):

t = np.array([m_i[i]])

cross_p = (t.T).dot(t)

A += cross_p * n_i[i]

# The within-group scatter matrix

W = np.zeros((n_dim, n_dim))

for i in range(n_groups):

t = Z[y == i + 1] - m_i[i]

for v in t:

vr = np.array([v])

W += (vr.T).dot(vr)

from eigen_sym import eigen_sym as eigen

Lmbd, U = eigen(W) # Eigen vectors in row format

# U[i]: the i-th eigen vector

L_sqrt = Lmbd ** 0.5

ULsqrt = (U.T) @ np.diag(L_sqrt)

ULsqrt_inv = np.linalg.inv(ULsqrt)

ULAUL = ULsqrt_inv @ A @ (ULsqrt_inv.T)

LmbdA, ULv = eigen(ULAUL)

print('Lambda =')

fout.write('\nLambda = \n')

for tv in LmbdA:

print(' ', tv)

fout.write(' {}\n'.format(tv))

plt.bar(range(n_dim), LmbdA, align = 'center')

plt.xlabel('$\lambda$', fontsize = 14)

plt.ylabel('Eigen value')

plt.show()

from scipy.stats import chi2

for k in range(n_groups - 1):

print('\nk = ',k)

fout.write('\nk = {}\n'.format(k))

df = (n_dim - k) * (n_groups - k - 1)

lambda_p = 1.0

for j in range(k, n_groups):

lambda_p *= 1.0 + LmbdA[j]

v_chi2 = -(n - (n_dim + n_groups)/2 - 1) * np.log(1 / lambda_p)

cum_chi2 = chi2(df).cdf

p_value = 1.0 - cum_chi2(v_chi2)

print('p = {0:<9.5f} chi2 = {1:<15.3f} df = {2:}'.format(p_value, v_chi2, df))

fout.write('p = {0:<9.5f} chi2 = {1:<15.3f} df = {2:}\n'.format(p_value, v_chi2, df))

v = np.linalg.inv(ULsqrt.T) @ (ULv.T)

v = v.T

for i in range(len(v)):

tv = np.array([v[i]])

v_sum = (tv @ W @ tv.T)[0][0]

v[i] = v[i] / (v_sum ** 0.5)

print('\nv =')

fout.write('\n\nv =\n')

for vv in v[:n_groups - 1].T:

for i in range(len(vv)):

print('{0:>9.5f} '.format(vv[i]), end = '')

fout.write('{0:>9.5f} '.format(vv[i]))

print()

fout.write('\n')

yj = np.empty((n_groups - 1, n))

for j in range(n_groups - 1):

yj[j] = np.array([v[j]]) @ Z.T

T = np.zeros((n_dim, n_dim))

for vv in Z:

ndv = np.array([vv])

T += (ndv.T).dot(ndv)

thetas = np.empty(n_groups - 1)

for j in range(n_groups - 1):

thetas[j] = (np.array([v[j]]) @ (T * (1/n)) @ (np.array([v[j]]).T))[0][0]

fj = np.empty((n_groups - 1, n))

for j in range(n_groups - 1):

fj[j] = yj[j] * (1 / (thetas[j]**0.5))

R = np.empty((n_dim, n_groups - 1))

for j in range(n_dim):

for g in range(n_groups - 1):

R[j][g] = ss.pearsonr(Z[:,j].T, fj[g])[0]

print('\nStructure matrix')

fout.write('\n\nStructure matrix\n')

print('{0:<10s}'.format('Variable'), end = '')

fout.write('{0:<10s}'.format('Variable'))

for j in range(n_groups - 1):

print('{0:>9s}'.format('Func-{}'.format(j + 1)), end = '')

fout.write('{0:>9s}'.format('Func-{}'.format(j + 1)))

print()

fout.write('\n')

for i in range(n_dim):

print(' {0:<9s}'.format(VarID[i]), end = '')

fout.write(' {0:<9s}'.format(VarID[i]))

for j in range(n_groups - 1):

print('{0:>9.5f}'.format(R[i][j]), end = '')

fout.write('{0:>9.5f}'.format(R[i][j]))

print()

fout.write('\n')

dim1 = 1

dim2 = 2 if n_groups > 2 else dim1

while True:

plt.figure(figsize = (6, 6))

plt.plot(R.T[dim1 - 1], R.T[dim2 - 1], 'o')

pos_id = 0

for tx, ty in zip(R.T[dim1 - 1], R.T[dim2 - 1]):

plt.text(tx + 0.01, ty + 0.01, VarID[pos_id])

pos_id += 1

plt.xlabel('Function {}'.format(dim1))

plt.ylabel('Function {}'.format(dim2))

plt.xlim([-1.1, 1.1])

plt.ylim([-1.1, 1.1])

plt.plot([-1, 1], [0, 0], color = 'k')

plt.plot([0, 0], [-1, 1], color = 'k')

plt.title('Map of\nCorrelations of Variables and Discriminant Functions')

plt.tight_layout()

plt.show()

print('\nChoose discriminant functions for the Structure Map: from {0} to {1}'.

format(1, n_groups - 1))

dim1 = int(input('Dim1 = '))

if dim1 < 1:

break

if dim1 > n_groups - 1:

dim1 = n_groups - 1

dim2 = int(input('Dim2 = '))

if dim2 < 1:

dim2 = 1

if dim2 > n_groups - 1:

dim2 = n_groups - 1

markers = ['1', '2', '3', '4', '8', 's', 'p', '*', 'h', '+', 'x', 'D']

dim1 = 1

dim2 = 2 if n_groups - 1 > 1 else 1

while True:

if dim1 == dim2:

for g in range(n_groups):

sb.kdeplot(fj[dim1 - 1][y == g + 1], label = g + 1)

plt.xlabel('Function {}'.format(dim1))

plt.show()

else:

for g in range(n_groups):

plt.scatter(fj[dim1 - 1][y == g + 1], fj[dim2 - 1][y == g + 1],

marker = markers[g], label = (g+1))

plt.xlabel('Function {}'.format(dim1))

plt.ylabel('Function {}'.format(dim2))

plt.legend()

plt.show()

print('\nChoose discriminant functions for the Distribution Map...1 to {}'.

format(n_groups - 1))

dim1 = int(input('Dim 1 = '))

if dim1 < 1:

break

elif dim1 >= n_groups:

dim1 = n_groups - 1

dim2 = int(input('Dim 2 = '))

if dim2 < 1:

break

elif dim2 >= n_groups:

dim2 = n_groups - 1

fout.close()

print('\n\nOutput file <{}> was saved.\n'.format(fout_name))

上のプログラムを実行すると、入力データファイル名と出力用ファイル名の設定が求められる（図１）。

図１

入力データファイルは、CSV形式のファイルとして用意する。CSV形式のファイルは、Excelで簡単に用意できる（図２）。

図２

一行目に、変数の名前を設定し、２行目から１行に１ケースのデータを設定する。第１欄は、各ケースのデータを区別するための文字列を設定する。このケースを区別する文字列は、文字列ではなく数値を設定するとプログラムの実行時にエラーとなる。文字列に数字が含まれるのは問題ない。第２欄は、各ケースのカテゴリを整数値で設定する。カテゴリー数（グループ数）がｇであれば、カテゴリーは１からｇまでの整数値で表す。第３欄から各ケースの変数値を設定する。図２では、X1からX5までの５個の変数の値が設定されている。Excelでのデータの設定が終わったら、「名前を付けて保存」を選ぶと、図３のダイアログボックスが表示される。

図３

「ファイルの種類」として「CSV（カンマ区切り）（*.csv）」を選ぶと、CSV形式で保存される。

出力用ファイル名は適当なテキストファイル名を設定すればよい。出力用ファイル名を設定してEnterキーを押すと、計算が始まり、図４のグラフが表示される。

図４

図４は、固有値の大きさを示すものである。固有値は、０でないものはグループ数ｇより１少ない個である。上の例の場合、４グループなので、個の固有値が０でないことが示されている。描かれたグラフは、保存ボタンをクリックするとファイルに保存される。

図４のウィンドウの右上のX印をクリックしてウィンドウを閉じると、次に図５の散布図が表示される。

図５

図５は、変数と判別関数の相関係数を座標値として描かれたものである。横軸にFunction 1（判別関数 1）と変数との相関係数、縦軸にFunction 2（判別関数 2）と変数との相関係数を座標値とするものである。変数X1と判別関数１との相関係数が１に近く、変数X2と判別関数２との相関係数が－0.8ぐらいであることがわかる。判別関数は、求められた個の固有値に対応して判別関数 1から判別関数まで求められる。図５のフォームを閉じると、次に描く散布図の判別関数の選択が求められる（図６）。