アセットマネージャーのためのファイナンス機械学習：特徴量の重要度分析 ロジスティック回帰のp値とMDI /MDA
ロジスティック回帰とは、幾つかの説明変数から、２値からなる結果が起こる確率を説明、または予測する統計手段である。
説明変数を$${{\bf X},(X_1,\dots, X_n) }$$、結果を$${y \in \{0,1\}}$$とする。$${y=1}$$となる確率を$${P(y=1)}$$とすれば、$${y=0}$$となる確率は、$${1-P(y=1)}$$である。
ロジスティック回帰は、結果の$${y}$$をロジット変換し、説明変数$${{\bf X}}$$を線形関係に結びつけ、その係数$${{\bf \beta}, \beta_0, \dots \beta_n}$$を推定する。
二値変数$${y}$$の確率$${P(y=1)}$$をロジット変換$${\displaystyle{\log(\frac{P(y=1)}{1-P(y=1)})}}$$することにより、
$${\begin{cases} P(y=1)\rightarrow 0 & \displaystyle{\log(\frac{P(y=1)}{1-P(y=1)})}\rightarrow -\infty \\ P(y=1)\rightarrow 1 & \displaystyle{\log(\frac{P(y=1)}{1-P(y=1)})}\rightarrow +\infty \end{cases} }$$
の連続関数にして、
$${\displaystyle{\log(\frac{P(y=1)}{1-P(y=1)}) = \beta_0 + \sum^{n}_{i=1}\beta_i X_i}}$$
とできる。
一般に、$${\displaystyle{\log(\frac{P(y=1)}{1-P(y=1)})}}$$を$${y=1}$$が起こる確率のオッズと呼ぶ。
またロジット変換の右辺を$${t}$$と置き、逆変換をすれば、
$${t=logit(P(y=1))=\displaystyle{\log(\frac{P(y=1)}{1-P(y=1)}) }}$$から、
$${P(y=1)=\displaystyle{\frac{e^t}{1+e^t} = \frac{1}{1+e^{-t}}}}$$
が得られる。
これから説明変数が結果$${y=1}$$にどれくらいの影響を与えているかの偏回帰係数の推定は、$${y=1}$$が最も起こりやすくなるように観測値$${{\bf x}}$$にかかる係数$${{\bf beta}}$$を求める最尤法を使って求める。
２値の結果$${y}$$が$${m}$$個あり、互いに独立であるとすれば、$${y_i}$$が$${1}$$を取る確率を$${P(y_i=1)=P_i}$$と書けば、、$${P_i^{y_i}(1-P_i)^{1-y_i}}$$で与えられるから、$${{\bf y},(y_1, \dots, y_m)}$$の尤度は、
$${L=\Pi^{m}_{i=1}P_i^{y_i}(1-P_i)^{1-y_i}}$$
と書ける。
この対数を取り、
$${\log(L)=\displaystyle{\sum^{n}_{i=1}(y_i \log P_i + (1-y_i)\log (1-P_i)}}$$、
これを、対数尤度関数と呼ぶ。
ここに、$${P_i=\displaystyle{ \frac{1}{1+e^{-t}} }}$$ $${\displaystyle{=\frac{1}{1+e^{-(\beta_0 + \beta_1 x_1 + \dots + \beta_n x_n)}} }}$$を代入し、$${\log(L)}$$が最大になる最尤推定値$${{\bf \hat{\beta}}}$$を計算する。
偏回帰係数の$${\beta_i}$$について、$${\beta_i > 0}$$である場合は、説明変数$${x_i}$$は事象$${y=1}$$に対し、ポジティブ（事象が起こる方）に作用すると解釈でき、$${\beta_i < 0}$$である場合は、説明変数$${x_i}$$は事象$${y=1}$$に対し、ネガティブ（事象が起こらない方）に作用すると解釈できる。
また、$${\beta_i}$$についてのp値が0.05未満だった場合、有意性があり、$${\beta_i}$$がゼロである帰無仮説は棄却される。
このロジスティック回帰はstatsmodelに実装されている。
これを用いて、情報的な説明変数を５、冗長な説明変数を30、残りはノイズの40の説明変数$${{\bf X}}$$、結果$${{\bf y}}$$を$${0,1}$$として、サンプル数1000についてロジット回帰を行う。
import numpy as np
import pandas as pd
def getTestData(n_features=100, n_informative=25, n_redundant=25, n_samples=10000, random_state=0, sigmaStd=.0):
import FeatImp as fi
importlib.reload(fi)
#generate a random dataset for classification problem
from sklearn.datasets import make_classification
np.random.seed(random_state)
kwds=dict(n_samples=n_samples,n_features=n_features-n_redundant,
n_informative=n_informative,n_redundant=0,
random_state=random_state)
X, y = fi.getTestData(**kwds)
cols = ['I_'+str(i) for i in range(0, n_informative)]
cols += ['N_'+str(i) for i in range(0, n_features - n_informative - n_redundant)]
X = pd.DataFrame(X, columns=cols)
i = np.random.choice(range(0, n_informative), size=n_redundant)
for k, j in enumerate(i):
X['R_'+str(k)] = X['I_' + str(j)] + np.random.normal(size=X.shape[0])*sigmaStd
return X, y ここでは、ファイナンス機械学習の特徴量の重要度で使用したgetTestDataを呼び出している。
import statsmodels.discrete.discrete_model as sm
import seaborn as sns
import matplotlib.pylab as plt
X, y = getTestData(40, 5, 30, 10000, sigmaStd=.1)
ols=sm.Logit(y['bin'],X).fit()
plot_data = ols.pvalues.sort_values(ascending=True)
plot_data.plot(kind='barh', figsize=(20,10),
title="p-Values computed on a set of explanatory variables")
plt.show()
同じデータセットを、決定木のバギングでMDI分析する。
import FeatImp as fi
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import BaggingClassifier
clf = DecisionTreeClassifier(criterion='entropy',
max_features=1,
class_weight='balanced',
min_weight_fraction_leaf=0)
clf = BaggingClassifier(estimator=clf,
n_estimators=1000,
max_features=1.,
max_samples=1.,
oob_score=False)
fit = clf.fit(X,y['bin'])
imp = fi.getFeatImpMDI(fit, featNames=X.columns)
imp.sort_values('mean', inplace=True)
plt.figure(figsize=(10, imp.shape[0] / 5))
imp['mean'].plot(kind='barh', color='b', alpha=0.25, xerr=imp['std'], error_kw={'ecolor': 'r'})
plt.title(' MDI results')
plt.show()
同じデータセットに、同じ決定木バギングでMDAを適用する。
この記事が気に入ったらサポートをしてみませんか？