第7章 線形モデル編: 第一節 線形回帰モデル
こんにちは、やっと第7章まで来ましたね(*'ω'*)
このチャプターでは、様々な線形モデルについて取り扱います。この記事ではまず線形回帰モデルからスタートしたいと思います。
線形回帰ー入門
線形回帰は、連続的な目的変数と一つないし複数の説明変数（特徴量や独立な確率変数）を関係づけるもので、それらの関係を線形であると仮定します。すなわち以下を仮定します。
・一つの特徴量を取ってきたときに、その特徴量と目的変数の関係は直線である。
・この線の傾きは、他の変数の値には依存しません。
・各変数の出力は加法的である。（ただし、2つの変数の相互作用を表す新しい変数を含めることができます）
言い換えると、モデルは、この線形関係からのランダムな偏差を除いて、目的変数が特徴の線形結合によって説明または予測できることを前提としています。
インポートと設定
import warnings
warnings.filterwarnings('ignore')%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
from sklearn.linear_model import SGDRegressor
from sklearn.preprocessing import StandardScalersns.set_style('whitegrid')
pd.options.display.float_format = '{:,.2f}'.format単純回帰
まずは、ランダムなデータを生成します。ただし、少しだけ規則を持ちます。
x = np.linspace(-5, 50, 100)
y = 50 + 2 * x + np.random.normal(0, 20, size=len(x))
data = pd.DataFrame({'X': x, 'Y': y})
ax = data.plot.scatter(x='X', y='Y', figsize=(14, 6))
sns.despine()
plt.tight_layout()
左側に単一の独立変数を持つ線形モデルは、次のモデルを想定しています。
データが実際に直線に適合しない場合に発生する偏差またはエラーを考慮します。 𝜖はエラーまたは残差と呼びます。
statsmodelsを利用して単純回帰モデルを推定
サマリーの上部には、データセットの特性、つまり推定方法、観測値とパラメーターの数が表示され、標準誤差の推定では不均一性が考慮されていないことが示されます。
中央のパネルは、人工データ生成プロセスを厳密に反映する係数値を示しています。要約結果の中央に表示される推定値は、以前に導出されたOLS式を使用して取得できることを確認できます。
X = sm.add_constant(data['X'])
model = sm.OLS(data['Y'], X).fit()
print(model.summary())
'''
OLS Regression Results
==============================================================================
Dep. Variable: Y R-squared: 0.729
Model: OLS Adj. R-squared: 0.726
Method: Least Squares F-statistic: 263.9
Date: Tue, 04 Aug 2020 Prob (F-statistic): 1.48e-29
Time: 05:05:36 Log-Likelihood: -434.79
No. Observations: 100 AIC: 873.6
Df Residuals: 98 BIC: 878.8
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 54.3613 3.256 16.695 0.000 47.900 60.823
X 1.9145 0.118 16.246 0.000 1.681 2.148
==============================================================================
Omnibus: 1.054 Durbin-Watson: 1.822
Prob(Omnibus): 0.590 Jarque-Bera (JB): 1.153
Skew: 0.208 Prob(JB): 0.562
Kurtosis: 2.677 Cond. No. 47.6
==============================================================================
'''ここで決定係数の確認を行います。こちらでは傾き(決定係数)が1.9145で切片が54.3613となっております。
beta = np.linalg.inv(X.T.dot(X)).dot(X.T.dot(y))
pd.Series(beta, index=X.columns)
'''
const 54.36
X 1.91
dtype: float64
'''一致しましたね(*'ω'*)
モデルと残差のプロット
data['y-hat'] = model.predict()
data['residuals'] = model.resid
ax = data.plot.scatter(x='X', y='Y', c='darkgrey', figsize=(14,6))
data.plot.line(x='X', y='y-hat', ax=ax);
for _, row in data.iterrows():
plt.plot((row.X, row.X), (row.Y, row['y-hat']), 'k-')
sns.despine()
plt.tight_layout();
重回帰モデル
説明変数が増えるだけです。以下2つの説明変数になります。
サンプルデータ生成
## Create data
size = 25
X_1, X_2 = np.meshgrid(np.linspace(-50, 50, size), np.linspace(-50, 50, size), indexing='ij')
data = pd.DataFrame({'X_1': X_1.ravel(), 'X_2': X_2.ravel()})
data['Y'] = 50 + data.X_1 + 3 * data.X_2 + np.random.normal(0, 50, size=size**2)
X = data[['X_1', 'X_2']]
y = data['Y']
## Plot
three_dee = plt.figure(figsize=(15, 5)).gca(projection='3d')
three_dee.scatter(data.X_1, data.X_2, data.Y, c='g')
sns.despine()
plt.tight_layout();
statsmodelsを利用して重回帰回帰モデルを推定
パネルの右上部分には、すべての係数がゼロで無関係であるという仮説を棄却するF検定とともに、先ほど説明した適合度の測定値が表示されます。同様に、t統計は、切片係数と両方の勾配係数が、当然のことながら非常に有意であることを示しています。
要約の下部には、残留診断が含まれています。左側のパネルには、正規性の仮説をテストするために使用されるスキューと尖度が表示されます。 OmnibusとJarque-Bera検定はどちらも、残差が正規分布しているという帰無仮説を棄却できません。Durbin-Watson統計は残差のシリアル相関をテストし、2に近い値を持ちます。2つのパラメーターと625の観測が与えられた場合、自己相関がないという仮説を棄却できません。
最後に、条件数は多重共線性に関する証拠を提供します。これは、入力データを含む設計行列の最大および最小固有値の平方根の比です。 30を超える値は、回帰に有意な多重共線性がある可能性があることを示しています。
X_ols = sm.add_constant(X)
model = sm.OLS(y, X_ols).fit()
print(model.summary())
'''
OLS Regression Results
==============================================================================
Dep. Variable: Y R-squared: 0.791
Model: OLS Adj. R-squared: 0.791
Method: Least Squares F-statistic: 1180.
Date: Tue, 04 Aug 2020 Prob (F-statistic): 2.01e-212
Time: 05:27:17 Log-Likelihood: -3327.4
No. Observations: 625 AIC: 6661.
Df Residuals: 622 BIC: 6674.
Df Model: 2
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 50.4612 1.991 25.350 0.000 46.552 54.370
X_1 1.0208 0.066 15.407 0.000 0.891 1.151
X_2 3.0522 0.066 46.070 0.000 2.922 3.182
==============================================================================
Omnibus: 1.886 Durbin-Watson: 2.021
Prob(Omnibus): 0.390 Jarque-Bera (JB): 1.698
Skew: 0.100 Prob(JB): 0.428
Kurtosis: 3.160 Cond. No. 30.0
==============================================================================
'''検算
beta = np.linalg.inv(X_ols.T.dot(X_ols)).dot(X_ols.T.dot(y))
pd.Series(beta, index=X_ols.columns)アウトプットを画像として書き込む方法
plt.rc('figure', figsize=(12, 7))
plt.text(0.01, 0.05, str(model.summary()), {'fontsize': 14}, fontproperties = 'monospace')
plt.axis('off')
plt.tight_layout()
plt.subplots_adjust(left=0.2, right=0.8, top=0.8, bottom=0.1)
モデルと残差のプロット
three_dee = plt.figure(figsize=(15, 5)).gca(projection='3d')
three_dee.scatter(data.X_1, data.X_2, data.Y, c='g')
data['y-hat'] = model.predict()
to_plot = data.set_index(['X_1', 'X_2']).unstack().loc[:, 'y-hat']
three_dee.plot_surface(X_1, X_2, to_plot.values, color='black', alpha=0.2, linewidth=1, antialiased=True)
for _, row in data.iterrows():
plt.plot((row.X_1, row.X_1), (row.X_2, row.X_2), (row.Y, row['y-hat']), 'k-');
three_dee.set_xlabel('$X_1$');three_dee.set_ylabel('$X_2$');three_dee.set_zlabel('$Y, \hat{Y}$')
sns.despine()
plt.tight_layout();
その他の検定方法はこちらを参照
確率的勾配降下法
確率的勾配下降法とは
確率的勾配降下法（かくりつてきこうばいこうかほう、英: stochastic gradient descent, SGD）とは、連続最適化問題に対する勾配法の乱択アルゴリズム。目的関数が、微分可能な和の形である事が必要。バッチ学習である最急降下法をオンライン学習に改良した物。
ここでの目的
SGDを用いて目的関数を定義し、上記で説明したOLS（最小二乗法）の結果が正しいのかを確認する。
ライブラリ
sklearnライブラリには、linear_modelsモジュールにSGDRegressorモデルが含まれています。この方法を使用して同じモデルのパラメーターを学習するには、勾配がスケールに敏感であるため、最初にデータを標準化する必要があります。
データの用意
SGDRegressorは単位に敏感であるので、正規化しておく。
scaler = StandardScaler()
X_ = scaler.fit_transform(X)パラメーターの初期化
sgd = SGDRegressor(loss='squared_loss',
fit_intercept=True,
shuffle=True,
random_state=42,
learning_rate='invscaling',
eta0=0.01,
power_t=0.25)フィッティング
sgd.fit(X=X_, y=y)
'''
SGDRegressor(random_state=42)
'''結果の確認
coeffs = (sgd.coef_ * scaler.scale_) + scaler.mean_
pd.Series(coeffs, index=X.columns)
'''
X_1 1,037.02
X_2 2,824.40
dtype: float64
'''resids = pd.DataFrame({'sgd': y - sgd.predict(X_),
'ols': y - model.predict(sm.add_constant(X))})resids.pow(2).sum().div(len(y)).pow(.5)
'''
sgd 46.37
ols 46.37
dtype: float64
'''ここでSGDと最小二乗法の結果が一致していることが確認出来ました。
resids.plot.scatter(x='sgd', y='ols')
sns.despine()
plt.tight_layout();
まとめ
今回は線形回帰モデルについて紹介しました。
次に、線形回帰を使用してマルチアセットプライシングモデルを推定します。
この記事が気に入ったらサポートをしてみませんか？