「データ指向のソフトウェア品質マネジメント」の「テスト工程での不具合検出数の予測」をPythonで行う
import math
import pandas as pd
from pandas import DataFrame, Series
import matplotlib.pyplot as plt
%matplotlib inline
まずは対象データを読み込む。これは書籍のサポートページからダウンロードできる。 ダウンロードに書籍記載のパスワードが必要なので生データをここには載せないが、 データ自体はFenton (2008)のものである。
ソースコード行数KLOCと不具合検出数Failures以外はカテゴリーデータ(順序尺度)で、D1は開発メンバーの経験、P5はステークホルダーの関与の度合い、P6は顧客とのコミュニケーション、P7は協力会社のマネジメントを表す。
data1 = pd.read_excel(u"4.2節/4-2-data.xls", index_col='Prj')
data1.columns
def scatter_corr_matrix(df):
pd.scatter_matrix(df)
plt.show()
return df.corr()
scatter_corr_matrix(data1)
import statsmodels.formula.api as smf
smf.ols(formula="Failures ~ KLOC", data=data1).fit().summary()
smf.ols(formula="np.log(Failures) ~ np.log(KLOC)", data=data1).fit().summary()
data1_1 = pd.read_excel(u"4.2節/4-2-data.xls", index_col='Prj', sheetname='data1-1')
data1_1.describe()
ソースコード行数がすでに説明変数に含まれているため、 不具合検出密度と各順序尺度の相関(スピアマンの順位相関係数)を計算する。
from scipy.stats import spearmanr
data1_1_d = data1_1.assign(Density = data1_1['Failures'] / data1_1['KLOC'])
data1_1_d.corr(method='spearman')
不具合検出密度と各順序尺度の関係を箱ひげ図で表す。
data1_1_d.boxplot(column='Density', by='P7', return_type='dict')
data1_1_d.boxplot(column='Density', by='D1', return_type='dict')
data1_1_d.boxplot(column='Density', by='P5', return_type='dict')
data1_1_d.boxplot(column='Density', by='P6', return_type='dict')
plt.show()
これを見ると不具合検出密度と相関が高いのはP7(協力会社のマネジメントが高いと不具合検出密度が低い)だが、「データ志向の~」ではP7はデータ件数が少ない(n=7)こと、ソースコード行数との相関係数も極めて大きいこと(多重共線性の疑い)を理由に退けている。
重回帰分析の説明変数をステップワイズで選択¶
P7以外の順序尺度をステップワイズで選択する。まずは全部を使う。
smf.ols(formula="np.log(Failures) ~ np.log(KLOC) + D1 + P5 + P6", data=data1).fit().summary()
ここから1つずつ減らしていく。
formulas = [
"np.log(Failures) ~ np.log(KLOC) + D1 + P5 + P6",
"np.log(Failures) ~ np.log(KLOC) + D1 + P5",
"np.log(Failures) ~ np.log(KLOC) + D1",
"np.log(Failures) ~ np.log(KLOC)",
]
steps = [[res.llf, res.aic, res.rsquared, res.rsquared_adj, formula]
for (res, formula) in [(smf.ols(formula=formula, data=data1).fit(), formula) for formula in formulas]]
DataFrame(steps, columns=['Log-Likelihood', 'AIC', 'R-squared', 'Adj. R-squared', 'formula'])
AICおよびAdj. R-squaredで判断するとlog(KLOC)とD1を使うモデルが最も良い。
smf.ols(formula="np.log(Failures) ~ np.log(KLOC) + D1", data=data1).fit().summary()
このサマリーを見るとD1のダミー変数の中でも2_Lと3_MのP値が高い。そこで「データ指向の~」では3_M以下の水準を統合することを提案している。
d1_map = {'1_VL': '3_M', '2_L': '3_M', '3_M': '3_M', '4_H': '4_H', '5_VH': '5_VH'}
data1_d1_new = data1.assign(D1_new = data1['D1'].apply(lambda x: d1_map[x]))
ols_final = smf.ols(formula="np.log(Failures) ~ np.log(KLOC) + D1_new", data=data1_d1_new).fit()
ols_final.summary()
scatter_corr_matrix(
pd.concat(
[data1.dropna(subset=['KLOC']).reset_index(drop=True),
Series(np.exp(ols_final.predict()), name='predict')],
axis=1)[['Failures', 'predict']])
上記では散布図も表示している。なお、「データ指向の~」(2015年1月13日第4刷)に掲載の散布図はX軸とY軸が逆のように思える(軸ラベル上Xが予測でYが実績なのだが、右上でなく左下の散布図になっている)。
相関係数が0.887であることから、うまく予測できているように見えると判断している。
このくだりにはちょっとよくわからないことろがあって、そもそも予測値と実績値の相関係数というのは予測モデルの評価の仕方として適当なのだろうかという疑問がある。 というのはモデルが誤っていて「本来の値の約半分の値を予測する」というものであった場合、予測モデルとしては役に立たないが、相関係数は半分でない場合と同様に高くなるはずだからだ。
次に予測誤差 MRE (Megnitude of Relative Error)と呼ばれる評価指標を出している。こちらのほうが直感的な定義に思える。MREは1つ1つの予測に対する値なので、全体で1つの指標値にするには予測誤差の平均とか中央値とかで見るようだ。
def mre(actual, predicted):
df = DataFrame({'actual': actual, 'predicted': predicted})
mre = df.eval("abs(actual - predicted) / actual")
mre.plot.box()
plt.show()
print mre.describe().round(2)
return mre.median()
mre(
data1.dropna(subset=['KLOC'])['Failures'],
np.exp(ols_final.predict()))
import statsmodels.api as sm
poisson_model = smf.glm(
formula='Failures ~ np.log(KLOC) + D1 + P5 + P6',
data=data1,
family=sm.families.Poisson()).fit()
poisson_model.summary()
scatter_corr_matrix(
pd.concat(
[data1.dropna(subset=['KLOC']).reset_index(drop=True),
Series(poisson_model.predict(), name='predict')],
axis=1)[['Failures', 'predict']])
予測値と実測値の相関係数は向上している。
一方MREを計算すると中央値では改善するが平均値では悪くなっており、またばらつきも大きくなっている。
mre(
data1.dropna(subset=['KLOC'])['Failures'],
poisson_model.predict())