Up

官能評価に使える統計解析ツール Python

日本官能評価学会誌、第26巻第2号、102-106.

 

 

 本コラムでは、Pythonの使い易さを幾つか例示したが、その1つとして特異値分解(Singular Value Decomposition:SVD)を取り上げ、CATAデータのSVDによる分析例を示した。

 特異値分解SVDでは、行列を次式で展開する(岡本、2019p. 67岡本、2009p.47)。

データサイエンスでは、右辺の展開項の一部を用いて情報の圧縮が行われる(Brunton & Kutz, 2019)。ベクトルは、行列の行に関する情報と列に関する情報を表すものであり、は重みである。SVDCATA(Check-All-That-Apply method)データの分割表(クロス表)に適用してみる。各ベクトルの平方根で重み付けたものを座標値として図示すると行および列に含まれる関係の情報が視覚的に読み取れる。行の座標値と列の座標値を重ねて図示すると行と列の関係が1つの図に表され、関係を読み取るのに便利である。SVDを内積として解釈することは、追記で説明した。内積として解釈すると、クラスター分析を行うことができるが、クラスター分析の結果も示した。

CATAデータとして島村ら(2017Table 2を用いてみる(Excelに設定したものを図1に示す)。

図1

 

図1のデータ(ExcelファイルCATAData.xlsx)をSVDで分析するスクリプトCATADataSVD.pyリスト1のように作成した。スクリプトファイルCATADataSVD.pyとサンプルデータファイルCATAData.xlsxは圧縮ファイルsamplefiles.zipにまとめた。ダウンロード・展開して自由に利用されたい。

リスト1のスクリプトを実行すると、図2のグラフが描かれる。

図2

 

第1軸(横軸)の値で見ると、米は、コシヒカリ、ミルキークイーン、ササニシキのグループと低たんぱく米とに分かれている。評価語は、「好きな味」と「程良い粘り」がまとまっており、コシヒカリ、ミルキークイーン、ササニシキのグループと近い位置である。これは、図1における「好きな味」と「程良い粘り」の頻度がコシヒカリ、ミルキークイーン、ササニシキにおいて高いことに対応している。第2軸(縦軸)の座標について見てみると、値が0から離れているものとしては、米ではササニシキが正の側、低たんぱく米が負の側にある。評価語を見ると、正の側に「好きな味」があり、ササニシキに近い位置である。負の側には、「嫌いな味」、「ねっとりした」、「渋味」があり、低たんぱく米の近くに位置している。図3におけるこれらの組み合わせに対する分割表の頻度の値が高いことに対応している。

 なお、図2において横軸の値が全て負であるが、特異値分解の場合、ベクトルが特異ベクトルの組であれば、それらの符号を同時に反転したベクトルも特異ベクトルの組であるので、第1軸において、米と評価語の値を同時に反転すれば、第1軸の値が正である布置を得る。

 

追記

岡太彬訓 立教大学名誉教授より、式(1)を行列の成分の内積による表現として見るという観点のコメントを頂いた(私信、2022.11.02)。この観点では、式(1)を以下のように展開する。

特異値分解を次式で表す。

ここで、

, ,

である。

いま、

, 

と置くと、

である。すなわち、

を得る。

クロス表の各セルの値が内積で与えられる。これらの内積を与えるベクトルは、同じ空間内のベクトルである。所与のデータは、内積であると解釈されるモデルである。ベクトル間の距離は、そのモデルから導かれるもので、距離に対応するデータはない。

式(2)を変形して、

, 

を得る。内積モデルにおいては、座標の直交回転ができる。

内積モデルで考えると、クラスター分析が可能である。

 

クラスター分析

内積を与える座標ベクトルは、同じ空間内のベクトルであるので、これに基づいてクラスター分析を行うことができる。リスト1のスクリプトでは、クラスター分析も行っている。リスト1のスクリプトを実行すると、図3の結果が表示される。

図3

 

なお、統計分析のPythonプログラム例は、<岡本安晴「データを使いこなすための統計入門」22世紀アート、電子書籍>のホームページにも用意している。

 

 

引用・参考文献

Brunton, S. L., & Kutz, J. N. (2019). Data-driven science and engineering. New York: Cambridge University Press.

岡太彬訓(私信、2022.11.02

岡本安晴(2009) 統計学を学ぶための数学入門[下].培風館

岡本安晴(2019) いまさら聞けないPythonでデータ分析.丸善出版

島村 綾、小泉昌子、峯木眞知子、市原 茂(2017) 飯の官能評価の時系列変化.日本家政学会誌、689478-485

 

 

リスト1 CATAデータのSVDCATADataSVD.py

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

from scipy.cluster.hierarchy import linkage

from scipy.cluster.hierarchy import dendrogram

 

#   ExcelファイルをPandasDataFrame型として読み込む

xlfile = pd.ExcelFile('CATAData.xlsx')

df = pd.read_excel(xlfile)

 

#    1行目の値(米の種類)

col_names = df.keys()         #    Excelの1行目

col_names = col_names[1:]     #    2番目以降を取り出す

 

#    '評価用語'の列を取り出す

row_names = df['評価用語']

 

#    セル内の数値

df_data = df.values           #   2行目以降が取り出される

df_data = df_data[:,1:]       #    取り出した各行の2列目以降を取り出す

 

#    セル内の数値を実数型に変換

data = np.asfarray(df_data)

 

#     dataの特異値分解SVD

u, s, vh = np.linalg.svd(data, full_matrices = False)

 

print('\ns =', s)

 

#     特異値の平方根による重み付け

col_coord = np.diag(s**0.5) @ vh

row_coord = u @ np.diag(s**0.5)

 

#     日本語表示用

plt.rcParams['font.family'] = ['Yu Gothic',

'Noto Sans CJK JP','Hiragino Maru Gothic Pro']

 

#

#     重み付けられた特異ベクトルによる描画

#

plt.plot(row_coord.T[0],row_coord.T[1],'^',c = 'k')

plt.plot([], '^',c = 'k',label = '評価項目')

plt.plot(col_coord[0],col_coord[1],'v',c = 'k')

plt.plot([], 'v', c = 'k', label = '')

for i in range(len(col_names)):

    plt.text(col_coord[0][i], col_coord[1][i],

             col_names[i])

for i in range(len(row_names)):

    plt.text(row_coord.T[0][i], row_coord.T[1][i],

             row_names[i],rotation = -20,

             rotation_mode = 'anchor')

plt.legend()

plt.savefig('FigMap.png')

plt.show()

 

#

#       クラスター分析のためにクロス表の行と列の座標を統合

#

mrg_coord = np.concatenate((row_coord, col_coord.T), axis = 0)

mrg_label = np.concatenate((row_names, col_names), axis = 0)

#

#       クラスター分析

#

row_clusters = linkage(mrg_coord, method = 'ward')

fig = plt.figure(figsize = (15, 8), facecolor = 'white')

row_dendr = dendrogram(row_clusters, labels = mrg_label,

                       orientation = 'right')

plt.savefig('FigCluster.png')

plt.show()

 

 

Up