Home

拙著「Pythonでデータ分析」の第８章「数量化」の更新と補足

数量化のスクリプトをExcelファイルを読み込む方法に変更するとともに、スクリプトの更新を行った。変数の変換結果の出力もExcelファイルへの出力とした。

スクリプトはリスト１に示すとともに、ファイル類は圧縮ファイルc08files.zipにまとめたので、ダウンロード展開して自由に利用されたい。

外れ値を混合モデルによりベイズ分析を行う例を、本サイトの後半に示した。

入力データファイルのExcelファイル例を図１に示す。これは、拙著の図8.2.1のExcelファイル例と同じ内容である。拙著の図8.2.1ではCSVファイルとして保存する場合の説明が行われているが、ここでは、Excelファイルとして保存されるものとする。

図１

表の形式は、１行目が変数名、２行目が変換方法の指定で、０は無変換、１は標準化、２以上は当該変数のカテゴリ数で、カテゴリは１以上の整数値である。使われていないカテゴリがあれば、そのカテゴリは無視して、連続したカテゴリの整数値は全て使われているものとする。

リスト１のスクリプトを実行すると、入力データファイル名を聞いてくる。例えば、図１のファイルDataQ.xlsxを入力ファイルとするときは、次のようにファイル名を設定する。

データファイル(*.xlsx) = DataQ.xlsx

入力データファイル名を設定すると、計算が始まる。

計算が終了すると、以下のメッセージが表示される。

Results.xlsx was saved.

Results.txt was saved.

ファイルResults.xlsxには、変換後の値が保存されている（図２）。

図２

変数名の約束は、拙著p.136を参照されたい。

出力ファイルResults.txtを開くと、以下のような内容である。

Data file = DataQ.xlsx

Q1-0:

Weight =

0.40825

-0.81650

0.40825

Q1-1:

Weight =

0.70711

0.00000

-0.70711

Q2-0:

Weight =

-0.61920

0.77052

-0.15132

　　　・

　　　・

　　　・

Q5-1:

Weight =

-0.15960

-0.61367

0.77327

数量化における数量化ベクトルがWeightとして出力されている。拙著pp.137-138に記載されているものと符号が逆のものがあるが、これは固有ベクトルの符号を逆にしたものも解となるので、いずれも正しい。

数量化後の主成分分析の結果を図３に示す。

図３

数量化の変数値は、符号が反転したものも正しい値であることを考えると、図３は拙著の図8.2.4に対応している。

まず、赤枠で囲まれた変数に注目する。図３において、近い変数同志は同じ様に変化していること、逆に原点に関して反対側にある変数は逆に変化していることが示されている。拙著p.140の図8.2.5（Q6とQ2-1）および図8.2.6（Q6とQ4-0）と同じく単回帰分析を行うと次図の結果を得た。

図４

拙著の図8.2.5とは横軸が逆方向になっているが、同じ布置である。変数値が逆符号のものも正しい値である。

図５

図５は、拙著図8.2.6と同じである。

図３の布置を見るとQ1-1とQ2-0が重なっている（拙著図8.2.4では変数の符号が逆であるので、原点を挟んで反対側に位置している）。これら2つの項目の回答パターンの関係を調べるために、クロス表を作成してみる。クロス表は、次のスクリプトで簡単に作成できる。

import pandas as pd

xlsx = pd.ExcelFile('DataQ.xlsx')

data_xlsx = pd.read_excel(xlsx)

xtable = pd.crosstab(data_xlsx['Q1'][1:], data_xlsx['Q2'][1:])

print(xtable)

実行すると以下の出力を得る。

Q2 1 2 3

Q1

1 1 10 1

2 3 10 2

3 6 4 2

Q1-1の数量化ベクトルは

（0.7、 0、 -0.7）

であり、Q2-0の数量化ベクトルは

（-0.6、 0.8、 -0.2）

である。数量化ベクトルの成分の大きさの順に並べ替えるため、上のクロス表の行の並べ替えのインデックスのリストを

［２，１，０］

列の並べ替えのインデックスを

［０，２，１］

と置き、行と列の並べ替えを次のスクリプトで行う（並べ替えについてはこのウェブサイトで説明している）。

v_table = xtable.values

print(v_table)

print(v_table[[2,1,0]])

print(v_table[:,[0,2,1]])

実行すると以下の出力を得る。

[[ 1 10 1]

[ 3 10 2]

[ 6 4 2]]

[[ 6 4 2]

[ 3 10 2]

[ 1 10 1]]

[[ 1 1 10]

[ 3 2 10]

[ 6 2 4]]

最後の出力は、拙著の表8.2.2と同じ内容である。

リスト１　数量化のスクリプト

#

# Common Quantification

# of

# Categorical Variables

#

# Yasuharu Okamoto, 2016.09, 2023.05

#

import pandas as pd

import math

import numpy as np

import scipy.stats as ss

in_fl_nm = input('データファイル(*.xlsx) = ')

with pd.ExcelFile(in_fl_nm) as f:

data = pd.read_excel(f)

f_out = open('Results.txt', 'w')

f_out.write('Data file = '+ in_fl_nm)

print(data)

var_names = [v for v in data.keys()]

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

var_cond = data.values[0] # 変換の種類

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

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

X = data.values[1:] # ２行目以降を取り出す

print('X =\n', X)

n_data, n_vars = np.shape(X)

print('n_data =', n_data, ' n_vars =', n_vars)

dict_v = dict()

for i_var in range(0, n_vars):

if (var_cond[i_var] < 0): # 不適当な値

print("var_cond < 0")

raise "var_cond error."

else:

if (var_cond[i_var] == 0): # 無変換

dict_v[var_names[i_var]] = X.T[i_var]

elif (var_cond[i_var] == 1): # 標準化

vX = ss.zscore(ToFloat(X.T[i_var]))

dict_v[var_names[i_var]] = vX

else: # 数量化

vX = X.T[i_var]

K = var_cond[i_var]

G = np.zeros((n_data, K))

for j in range(n_data):

G[j][int(vX[j])-1] = 1.0

f = np.ones((1,n_data)) @ G / n_data

H = G - np.ones((n_data, 1)) @ f

Lmbd, Wght = np.linalg.eigh(H.T @ H) # 拙著式(8.1.3)

idx = np.argsort(-Lmbd) # 大きいものから小さいものへの並べ換え用インデックス

Lmbd = Lmbd[idx][:-1]

Wght = Wght.T[idx][:-1].T

coeff = (Lmbd ** (-0.5)) * (n_data**0.5)

z = H @ Wght @ np.diag(coeff) # 拙著式(8.1.4)

for i in range(len(Lmbd)):

dict_v[f'{var_names[i_var]}-{i}'] = z.T[i]

f_out.write(f'\n{var_names[i_var]}-{i}:\n')

f_out.write(f' Weight =\n')

for j in range(K):

f_out.write(f' {Wght.T[i][j]:.5f}\n')

#

# 変換データのExcelファイルへの書き出し

#

with pd.ExcelWriter('Results.xlsx') as f:

pd.DataFrame(dict_v).to_excel(f, index=False)

print('\nResults.xlsx was saved.')

f_out.close()

print('Results.txt was saved.')

外れ値の混合モデルによるベイズ分析

上の図４あるいは図５に見られるように、データが回帰モデルに従うものと回帰モデルからの予測とは離れた分布に従うものがあるとき、混合モデル（mixture model; Martin, 2018）で表すことができる。すなわち、回帰モデルに従うデータと回帰モデルとは別の分布に従うデータの混合と考える。

回帰モデルを

で表し、回帰モデルとは離れた分布を

で表す。データが回帰モデルに従う確率をと置けば、回帰モデルではなく分布（B2）に従う確率はであるので、混合モデルを次式（B3）で設定する。

混合モデルのStanスクリプトをリストB1に示す。pystan 3.7.0対応である。なお、スクリプトファイルなどは、まとめて圧縮ファイルbysnfiles.zipとした。ダウンロード・展開して自由に利用されたい。

ベイズ分析のStanスクリプトは、Anaconda/Ubuntu/WSL/WindowsのPython 3.11仮想環境において実行した。その実行後、ＭＣＭＣサンプルの分析をAnaconda/WindowsのPython 3.11仮想環境において行った。

単一の単回帰モデルによる分析と比較するために、単回帰モデルによるベイズ分析のStanスクリプトをリストB2のように用意した。リストB2のStanスクリプトを用いてQ6とQ2-1の間の単回帰分析を行うスクリプトをリストB3、リストB4のように用意した。Q6とQ4-0の間の単回帰分析は、後で行う。リストB3の実行後、リストB4を実行すると、図B1のグラフが表示される（日本語部分は上から貼付けたものである）。

図B1

図B1のグラフは図４とほぼ同じである。

リストB１の混合モデルによるMCMCサンプリングを行うリストＢ５のスクリプトの実行後、ＭＣＭＣサンプルの分析を行うリストＢ６のスクリプトを実行すると図Ｂ２のグラフが表示される。

図B2

図B2において、緑の点はが0.5より大きい、したがって回帰モデルに従う傾向の方が強いもの、赤の点はが0.5以下のもの、したがって外れ値の分布の方に従う傾向の強いものを表している。

図B2を閉じると、図B3のグラフが表示される。

図B3

３８人のの事後分布を重ねて表したものである。０より小さい領域、および１より大きい領域に分布が広がっているのは、KDEによる密度関数であることによる。の事後分布は、０あるいは１に明瞭に分かれていることが分かる。

単回帰モデルと混合モデルの比較を行うため、リストB7のスクリプトにより係数b0の事後分布を描くと図B4のようになる。

図B4

混合モデルの方がやや左寄りであるが、ほとんど重なっている。

係数b1の事後分布は、図B5のようになっている。

図B5

単純な回帰モデルの方が混合モデルより負の方向に分布がずれていることが見られる。「外れ値」を考慮した回帰モデルの方が、独立変数の効果が小さいという傾向があるようである。

Q6とQ4-0の間の単回帰分析

Q6とQ4-0の間の単回帰分析も同様に行う。スクリプトはExcel入力ファイルが異なるだけであるのでサイト上には示さないが、圧縮ファイルbysnfiles.zipには含まれているので、ダウンロード展開して参照されたい。

単回帰モデルによる分析結果は図B6のようになる。

図B6

混合モデルによる分析結果は図Ｂ７のようになる。

図B7

の事後分布は図Ｂ８のようである。

図B8

分布が０あるいは１の付近に集中しており、回帰モデルに従うか従わないかがはっきり分かれている。

単回帰モデルと混合モデルの比較は、以下のようである。

係数b0の分布は、ほぼ重なっている（図B9）。

図B9

係数b1の分布（図B10）は、単回帰モデルの場合は負の領域に分布しているので、単回帰直線は独立変数の値が右に行くほど減少していると言える。しかし、混合モデルの場合は、単回帰モデルの係数より分布が右寄りにあるので、外れ値を考慮した場合、独立変数の値が右に行ったときの単回帰直線の減少はより緩やかであると言えよう。

図B10

参考文献

Martin, O. (2018). Bayesian analysis with Python, second edition. Packt

リストB1．混合モデルのStanスクリプト

functions {

real mix_p_lpdf(real y, real pi, real mu, real sgm, real mu1, real sgm1 ) {

real v;

v = pi*exp(normal_lpdf(y|mu, sgm)) +

(1.0 - pi)*exp(normal_lpdf(y|mu1, sgm1));

return log(v);

}

}

data {

int n;

vector[n] y;

vector[n] x;

}

parameters {

real b0;

real b1;

real<lower=0> sgm;

vector[n] thetai;

}

transformed parameters {

vector[n] pi;

real mu1;

real sgm1;

pi = inv_logit(thetai);

mu1 = -100.0;

sgm1 = sgm;

}

model {

thetai ~ normal(0.0, 10.0);

b0 ~ normal(0, 10);

b1 ~ normal(0, 10);

sgm ~ exponential(0.2);

for (i in 1:n) {

y[i] ~ mix_p(pi[i], b0+b1*x[i], sgm, mu1, sgm1);

}

}

リストB2、単回帰モデルのStanスクリプト

data {

int n;

vector[n] y;

vector[n] x;

}

parameters {

real b0;

real b1;

real<lower=0> sgm;

}

model {

b0 ~ normal(0, 100);

b1 ~ normal(0, 100);

sgm ~ exponential(0.01);

y ~ normal(b0 + b1*x, sgm);

}

リストB3．MCMCサンプリングを行うスクリプト

import pandas as pd

import numpy as np

import stan

import pickle

import arviz as az

"""

Excelファイルからの読み込み

"""

in_file_name = 'Q2.xlsx'

xlsx = pd.ExcelFile(in_file_name)

data_xlsx = pd.read_excel(xlsx)

y = data_xlsx.values.T[1]

print(y)

x = data_xlsx.values.T[2]

print(x)

n = len(y)

with open('simple_reg.stan', 'r') as f:

sm = stan.build(f.read(), data = {'n':n, 'y':y, 'x':x})

fit = sm.sample()

d_frm = fit.to_frame()

i_data = az.from_pystan(posterior = fit, posterior_model = sm)

with open('temp.pkl', 'wb') as f:

pickle.dump((d_frm, i_data), f)

print(az.summary(i_data))

リストB4．MCMCサンプルの分析を行うスクリプト

import pandas as pd

import numpy as np

import pickle

import arviz as az

import matplotlib.pyplot as plt

import japanize_matplotlib

"""

Excelファイルからの読み込み

"""

in_file_name = 'Q2.xlsx'

xlsx = pd.ExcelFile(in_file_name)

data_xlsx = pd.read_excel(xlsx)

y = data_xlsx.values.T[1]

print(y)

x = data_xlsx.values.T[2]

print(x)

n = len(y)

with open('temp.pkl', 'rb') as f:

d_frm, i_data = pickle.load(f)

min_x = np.min(x)

max_x = np.max(x)

b0 = np.median(d_frm['b0'])

b1 = np.median(d_frm['b1'])

plt.plot([min_x, max_x], [b0 + b1*min_x, b0 + b1*max_x])

plt.title(f'b0 = {b0:<.3f}, b1 = {b1:<.3f}', fontsize = 20)

plt.plot(x, y, 'o')

plt.show()

with open('SimpleQ2b0b1.pkl', 'wb') as f:

pickle.dump((d_frm['b0'], d_frm['b1']), f)

リストB5．混合モデルによるMCMCサンプリングを行うスクリプト

import pandas as pd

import numpy as np

import stan

import pickle

import arviz as az

"""

Excelファイルからの読み込み

"""

in_file_name = 'Q2.xlsx' # input('入力データファイル(*.xlsx) = ')

xlsx = pd.ExcelFile(in_file_name)

data_xlsx = pd.read_excel(xlsx)

y = data_xlsx.values.T[1]

print(y)

x = data_xlsx.values.T[2]

print(x)

n = len(y)

with open('mix_reg.stan', 'r') as f:

sm = stan.build(f.read(), data = {'n':n, 'y':y, 'x':x})

fit = sm.sample()

d_frm = fit.to_frame()

i_data = az.from_pystan(posterior = fit, posterior_model = sm)

with open('temp.pkl', 'wb') as f:

pickle.dump((d_frm, i_data), f)

print(az.summary(i_data))

print(d_frm.keys())

for i in range(n):

v = d_frm[f'pi.{i+1}']

print(i, ' ', np.median(v))

リストB6．リストB5のスクリプトによるサンプルの分析スクリプト

import pandas as pd

import numpy as np

import pickle

import arviz as az

import matplotlib.pyplot as plt

import seaborn as sb

import japanize_matplotlib

pd.set_option('display.max_rows', 1000)

pd.set_option('display.max_columns', 1000)

"""

Excelファイルからの読み込み

"""

in_file_name = 'Q2.xlsx'

xlsx = pd.ExcelFile(in_file_name)

data_xlsx = pd.read_excel(xlsx)

y = data_xlsx.values.T[1]

print(y)

x = data_xlsx.values.T[2]

print(x)

n = len(y)

with open('temp.pkl', 'rb') as f:

d_frm, i_data = pickle.load(f)

print(az.summary(i_data))

min_x = np.min(x)

max_x = np.max(x)

b0 = np.median(d_frm['b0'])

b1 = np.median(d_frm['b1'])

plt.plot([min_x, max_x], [b0 + b1*min_x, b0 + b1*max_x])

n = len(x)

pi = np.empty(n)

for i in range(n):

pi[i] = np.median(d_frm[f'pi.{i+1}'])

plt.title(f'b0 = {b0:.3f}, b1 = {b1:.3f}', fontsize= 16)

for i in range(n):

plt.plot(x[i], y[i], 'o', c = 'g' if pi[i] > 0.5 else 'r')

plt.plot([], 'o', c = 'g', label = 'pi>0.5')

plt.plot([], 'o', c = 'r', label = 'pi<=0.5')

plt.legend(loc = 'center', fontsize = 16)

plt.show()

az.plot_trace(i_data)

plt.show()

for i in range(n):

pi_i = d_frm[f'pi.{i+1}']

sb.kdeplot(pi_i)

plt.title('Posterior Distributions of pi', fontsize = 16)

plt.xlabel('pi', fontsize = 14)

plt.tight_layout()

plt.show()

with open('MixQ2b0b1.pkl', 'wb') as f:

pickle.dump((d_frm['b0'], d_frm['b1']), f)

リストB7．係数の事後分布の比較

import pickle

import numpy as np

import matplotlib.pyplot as plt

import seaborn as sb

with open('SimpleQ2b0b1.pkl', 'rb') as f:

Sb0s, Sb1s = pickle.load(f)

with open('MixQ2b0b1.pkl', 'rb') as f:

Mb0s, Mb1s = pickle.load(f)

sb.kdeplot(Sb0s, label = 'Simple')

sb.kdeplot(Mb0s, label = 'Mix')

plt.title('Posterior Distribution of b0\nQ6 on Q4', fontsize = 16)

plt.legend(fontsize = 16)

plt.show()

sb.kdeplot(Sb1s, label = 'Simple')

sb.kdeplot(Mb1s, label = 'Mix')

plt.title('Posterior Distribution of b1\nQ6 on Q4', fontsize = 16)

plt.legend(fontsize = 16)

plt.show()

Home