In [4]:
#-*- encoding: utf-8 -*-
'''
Ouyoutoukei HW1
'''
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import statsmodels.api as sm
np.set_printoptions(precision=3)
pd.set_option('display.precision', 4)
In [5]:
# csvをインポート
df = pd.read_csv( 'odakyu-mansion.csv' )
# 基礎統計量を表示
print(df.describe())
家の向きは東, 西, 南, 北のダミー変数（0または1。南東の場合、南と東の両方に1）に分解し変換して、
In [6]:
# サンプルサイズ
data_len = df.shape[0]
# 家の向きはdummyに
df['d_N'] = np.zeros(data_len, dtype=float)
df['d_E'] = np.zeros(data_len, dtype=float)
df['d_W'] = np.zeros(data_len, dtype=float)
df['d_S'] = np.zeros(data_len, dtype=float)
for i, row in df.iterrows():
for direction in ["N", "W", "S", "E"]:
if direction in str(row.muki):
df.loc[i, 'd_{0}'.format(direction)] = 1
# 先頭10件を表示
print(df.head(10))
欠損値は平均値で置き換える。
In [7]:
df = df.fillna(df.mean())
In [8]:
# 定数項も加える
X = sm.add_constant(df[['time', 'bus', 'walk', 'area',
'bal', 'kosuu', 'floor', 'tf', 'd_N', 'd_E', 'd_S', 'd_W', 'year']])
# 普通の最小二乗法
model = sm.OLS(df.price, X)
results = model.fit()
# 結果を表示
print(results.summary())
となる。
p値を見ると、kosuu, floorがほとんど無関係であるように見える。
kosuuには外れ値が1つある（kosuu=2080）ので、それを除いてみる。
In [9]:
print(df.loc[161])
df = df.drop(161)
再び最小二乗法を実行すると、
In [10]:
X = sm.add_constant(df[['time', 'bus', 'walk', 'area', 'bal',
'kosuu', 'floor', 'tf', 'd_N', 'd_E', 'd_S', 'd_W', 'year']])
model = sm.OLS(df.price, X)
results = model.fit()
print(results.summary())
やはりkosuu, floorのp値が大きいので、説明変数から除くと、
In [11]:
X = sm.add_constant(df[['time', 'bus', 'walk', 'area',
'bal', 'tf', 'd_N', 'd_E', 'd_S', 'd_W', 'year']])
model = sm.OLS(df.price, X)
results = model.fit()
print(results.summary())
となる。さらに、balと南向き以外の方角のダミー変数を説明変数から除く。
In [12]:
X = sm.add_constant(df[['time', 'bus', 'walk', 'area', 'tf', 'year', 'd_S']])
model = sm.OLS(df.price, X)
results = model.fit()
print(results.summary())
p値が大きい南向きのダミー変数d_E、築年数yearも説明変数から除く。
In [13]:
X = sm.add_constant(df[['time', 'bus', 'walk', 'area', 'tf']])
model = sm.OLS(df.price, X)
results = model.fit()
print(results.summary())
p値が大きいtfを取り除く
In [14]:
X = sm.add_constant(df[['time', 'bus', 'walk', 'area']])
model = sm.OLS(df.price, X)
results = model.fit()
print(results.summary())
修正済みR^2 = 0.783, F統計量のp値3.96e-59 を見ると、「新宿駅からの乗車時間」, 「バスの乗車時間」, 「徒歩時間」, 「部屋の広さ」の4つで十分に住宅価格を説明できていると考えられる。
最小二乗法1〜6と比較しても、AIC・BICは殆ど変わらないか、改善している。
あとは残差を検討して、誤差項に関する諸仮定が満たされているかをチェックする。
In [15]:
# 回帰に使った変数だけを抜き出す
new_df = df.loc[:, ['price', 'time', 'bus', 'walk', 'area']]
# 説明変数行列
exp_matrix = new_df.loc[:, ['time', 'bus', 'walk', 'area']]
# 回帰係数ベクトル
coefs = results.params
# 理論価格ベクトル
predicted = exp_matrix.dot(coefs[1:]) + coefs[0]
# 残差ベクトル
residuals = new_df.price - predicted
# 残差をplot
fig, ax = plt.subplots(figsize=(12, 8))
plt.plot(predicted, residuals, 'o', color='b', linewidth=1, label="residuals distribution")
plt.xlabel("predicted values")
plt.ylabel("residuals")
plt.show()
# 残差平均
print("residuals mean:", residuals.mean())
In [16]:
print(new_df.loc[12] )
new_df = new_df.drop(12)
X = sm.add_constant(new_df[['time', 'bus', 'walk', 'area']])
model = sm.OLS(new_df.price, X)
results = model.fit()
print(results.summary())
In [17]:
# 説明変数行列
exp_matrix = new_df.loc[:, ['time', 'bus', 'walk', 'area']]
# 回帰係数ベクトル
coefs = results.params
# 理論価格ベクトル
predicted = exp_matrix.dot(coefs[1:]) + coefs[0]
# 残差ベクトル
residuals = new_df.price - predicted
# 残差をplot
fig, ax = plt.subplots(figsize=(12, 8))
plt.plot(predicted, residuals, 'o', color='b', linewidth=1, label="residuals distribution")
plt.xlabel("predicted values")
plt.ylabel("residuals")
plt.show()
# 残差平均
print("residuals mean:", residuals.mean())
最小二乗法6の結果に比べ、ばらつきが均等になった。
In [18]:
# 残差をplot
fig = plt.figure(figsize=(18, 10))
ax1 = plt.subplot(2, 2, 1)
plt.plot(exp_matrix['time'], residuals, 'o', color='b', linewidth=1, label="residuals - time")
plt.xlabel("time")
plt.ylabel("residuals")
plt.legend()
ax2 = plt.subplot(2, 2, 2, sharey=ax1)
plt.plot(exp_matrix['bus'], residuals, 'o', color='b', linewidth=1, label="residuals - bus")
plt.xlabel("bus")
plt.ylabel("residuals")
plt.legend()
ax3 = plt.subplot(2, 2, 3, sharey=ax1)
plt.plot(exp_matrix['walk'], residuals, 'o', color='b', linewidth=1, label="residuals - walk")
plt.xlabel("walk")
plt.ylabel("residuals")
plt.legend()
ax4 = plt.subplot(2, 2, 4, sharey=ax1)
plt.plot(exp_matrix['area'], residuals, 'o', color='b', linewidth=1, label="residuals - area")
plt.xlabel("area")
plt.ylabel("residuals")
plt.legend()
plt.show()
どの説明変数と残差の間にも特徴的な相関関係は見られない: 仮定5は満たす。
area変数だけ残差のばらつき方が異なるので、何らかの対策をとったほうが良い可能性がある（すみません、わかりません）。
当てはまりの良い回帰モデルを作ることが出来たが、残差の性質、特に等分散性の仮定を置いて良いのかについては問題が残った。等分散性の仮定に問題がある場合は、重みを付けて最小二乗法を使う必要があるので、慎重に考える必要がある。
In [ ]: