gaiaskyの技術メモ

ベイズ、機械学習などのデータサイエンスがテーマ。言語はpythonで。

中古マンション価格における駅・路線の影響分析をpymcでやってみる。

pymc python ベイズ

世の中には非常に優秀な人がいるもので、中古マンション価格について、そのスペックだけでなく、駅・路線の効果を分析した方がいる。（すごい！）

http://smrmkt.hatenablog.jp/entry/2014/11/17/234538

原著者の方は幾つかのモデルを試行錯誤的に検討されており、最終的なモデルは以下のようになっている。

\displaystyle \begin{align*} 価格_i &= b_0+b_1 距離_i + b_2 築年_i + b_3 部屋数_i + b_4 床面積_i + 駅_j \\ 駅_j &= r_j + \sum^K_{k=0} \delta_{j,k} 路線_k \end{align*}

こちらもstan+Rで実装されている。これは色々なデータに応用できそうなので、勉強のためにpymc3で実装した。 MCMCの実行においては、サンプリング手法としてNUTSを使うとともに、高速な近似手法であるADVI（自動微分変分法）を適用し、 計算速度と結果を比較した。

結果から言えば、MCMC（NUTS）による推定では、原著者とほぼ同じ推定結果が得られたが、4時間程度掛かった。 一方、ADIVは推定結果は少し（？）異なるが、2分程度で済んだ。従って、試行錯誤中はADVIで行い、 ある程度モデルが固まったらMCMCでしっかり推定する、と言う使い分けをしようと思う。

1. pymc3によるモデル実装

駅効果の部分が難しかったが、いろいろ試した結果、以下のような実装とした。

with pm.Model() as model4:
    # 家のスペック（広さ、駅からの距離など）の効果
    beta = pm.Normal('beta', mu=0, sd=100**2, shape=5)
    # 路線効果
    s_rt = pm.Uniform('s_rt', 0, 100**2)
    rt = pm.Normal('rt', mu=0, sd=s_rt, shape=len(df.train.unique()))
    # 駅単体の効果
    s_rs_base = pm.Uniform('s_rs_base', 0, 100**2)
    rs_base = pm.Normal('rs_base', mu=0, sd=s_rs_base, shape=len(df.station.unique()))
    # 駅は、駅単体と路線の影響を受ける
    s_rs = pm.Uniform('s_rs', 0, 100**2)
    rs = pm.Normal('rs', mu=rs_base+pm.math.dot(rt,data['ST'].T), sd=s_rs, shape=len(df.station.unique()))
    # マンションの平均価格
    price = beta[0] + \
            beta[1]*data['Distance'] + \
            beta[2]*data['From'] + \
            beta[3]*data['Room'] + \
            beta[4]*data['Space'] + \
            rs[data['S']-1]
    # 誤差
    eps = pm.HalfCauchy('eps', beta=100)
    #
    y_pred = pm.Normal('y_pred', mu=price, sd=eps, observed=data['Y'])

グラフィカルモデルは以下の通り。

pm.model_to_graphviz(model4)

f:id:gaiasky:20180819202935p:plain

2. MCMC（NUTS）による推定

データの読み込み部分から記載すると、以下のようになる。 私のMacbook Air(2011年モデル!)では、4時間ほど掛かった。

# データは原著者が公開しているgithubから。
# https://github.com/smrmkt/reins/tree/master/regression
df = pd.read_csv('mantions.csv')
df.dropna(inplace=True)

st = pd.read_csv('station_train.csv', header=None)

data = {
    'Distance': df.distance.values,
    'From': df['from'].values,
    'Room': df.room.values,
    'Space': df.space.values,
    'ST': st.values,
    'S': df.station.values,
    'T': df.train.values,
    'Y': df.price.values
}

# モデルは上述のため省略。

# MCMCの実行
with model4:
    trace4 = pm.sample(10000, tune=2000, start=pm.find_MAP(method='Powell'), step=pm.NUTS())

logp = -19,899, ||grad|| = 3,581.3: : 5001it [00:10, 489.30it/s]
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [eps, rs, s_rs, rs_base, s_rs_base, rt, s_rt, beta]
Sampling 2 chains: 100%|██████████| 24000/24000 [3:53:32<00:00,  1.98draws/s]

3. ADIVによる推定

ADVIは非常に高速。数分で終了。

with model4:
    inference = pm.ADVI()
    approx = pm.fit(n=60000, method=inference)
    trace4_advi = approx.sample(5000)

Average Loss = 18,925: 100%|██████████| 60000/60000 [02:27<00:00, 407.68it/s]
Finished [100%]: Average Loss = 18,925

4. 両者の比較

計算時間は実行パラメータに依るところもあるが、MCMC：4時間、ADVI：2分で、圧倒的にADVIが早い。

MCMCによる推定結果は以下の通り（路線、駅は上位・下位それぞれ5位表示）。 大体、原著者の結果と同じと思う。

print_effectiveness(df, trace4, upper_lower_station_n=5)
mean sd mc_error hpd_2.5 hpd_97.5 n_eff Rhat
切片 -1657.99 21.13 1.66e-01 -1700.18 -1617.23 12662.80 1.0
駅からの距離 -4.78 0.10 6.65e-04 -4.98 -4.59 20982.55 1.0
築年 0.85 0.01 8.71e-05 0.83 0.87 13096.16 1.0
部屋数 -0.42 0.19 1.51e-03 -0.79 -0.06 16784.44 1.0
床面積 1.93 0.15 1.25e-03 1.65 2.23 15127.33 1.0
mean sd mc_error hpd_2.5 hpd_97.5 n_eff Rhat
横浜高速鉄道ＭＭ線 19.66 2.84 0.08 14.06 25.20 1079.64 1.0
東急田園都市線 10.25 2.02 0.07 6.41 14.33 667.58 1.0
東急東横線 9.33 1.81 0.06 5.80 12.89 748.61 1.0
京浜東北・根岸線 5.85 2.55 0.04 0.78 10.79 3999.36 1.0
東急目黒線 5.38 4.53 0.04 -3.62 14.21 10849.40 1.0
... ... ... ... ... ... ... ...
鶴見線 -1.55 3.23 0.05 -7.57 5.15 4981.11 1.0
京王電鉄相模原線 -1.88 3.15 0.06 -8.09 4.29 2647.75 1.0
小田急電鉄多摩線 -4.21 2.87 0.06 -9.85 1.37 2188.33 1.0
京浜急行電鉄大師線 -4.46 2.45 0.06 -9.13 0.52 1514.44 1.0
金沢シーサイドＬ -7.45 2.28 0.06 -11.95 -2.94 1293.02 1.0

23 rows × 7 columns

mean sd mc_error hpd_2.5 hpd_97.5 n_eff Rhat
みなとみらい 35.10 1.68 0.07 31.80 38.39 508.95 1.00
武蔵小杉 25.14 1.43 0.07 22.31 27.92 344.11 1.01
横浜 17.89 1.45 0.07 15.13 20.71 352.20 1.01
日本大通り 17.63 1.92 0.07 13.90 21.42 590.00 1.00
馬車道 17.57 2.24 0.07 13.27 22.05 807.19 1.00
... ... ... ... ... ... ... ...
産業道路 -8.70 3.83 0.07 -16.07 -1.02 2863.14 1.00
小島新田 -9.45 2.61 0.07 -14.62 -4.34 1173.79 1.00
海の公園柴口 -9.95 3.18 0.07 -16.04 -3.57 1940.90 1.00
八景島 -11.39 3.75 0.07 -18.58 -3.92 2534.88 1.00
南部市場 -12.09 3.76 0.07 -19.18 -4.49 2750.36 1.00

133 rows × 7 columns

ADVIの結果は以下。MCMCの推定結果と比較すると、ADVI目線で言うと、

  • 家スペック効果では、切片が大きくことなる。また、床面積の効果が少し大きめ。
  • 路線効果も順位の傾向は近いが、係数は大分小さい。家スペックの切片が小さいこととの兼ね合いと思われる。東海道本線が上位に来ている。
  • 駅効果も上位は傾向が近いが、下位は少し違うようである。また、係数はこちらも小さい。
print_effectiveness(df,trace4_advi, upper_lower_station_n=5)
mean sd mc_error hpd_2.5 hpd_97.5
切片 -0.39 2.27e-01 3.31e-03 -0.83 4.35e-02
駅からの距離 -5.52 8.87e-02 1.25e-03 -5.69 -5.35e+00
築年 0.02 1.15e-04 1.46e-06 0.02 1.97e-02
部屋数 -0.17 8.60e-02 1.16e-03 -0.34 -1.59e-03
床面積 3.78 5.86e-02 8.64e-04 3.66 3.89e+00
mean sd mc_error hpd_2.5 hpd_97.5
横浜高速鉄道ＭＭ線 2.85 3.17 0.04 -3.34 9.15
東急東横線 1.46 2.23 0.03 -2.77 5.96
東急田園都市線 1.23 2.40 0.03 -3.42 5.98
東海道本線 1.19 3.61 0.05 -6.12 7.93
東急目黒線 0.91 3.85 0.05 -6.36 8.56
... ... ... ... ... ...
京浜急行電鉄逗子線 -1.33 3.62 0.05 -8.57 5.52
横浜市ブルーライン -2.04 1.82 0.03 -5.71 1.37
小田急電鉄小田原線 -2.22 2.54 0.03 -7.02 3.01
金沢シーサイドＬ -2.98 2.74 0.04 -8.48 2.21
京浜急行電鉄本線 -3.17 1.73 0.03 -6.53 0.22

23 rows × 5 columns

mean sd mc_error hpd_2.5 hpd_97.5
みなとみらい 13.35 3.67 0.05 6.30 20.54
武蔵小杉 12.73 2.16 0.03 8.51 16.98
日本大通り 11.07 3.30 0.05 4.61 17.39
横浜 10.63 1.40 0.02 7.92 13.24
たまプラーザ 10.46 2.26 0.03 6.36 15.10
... ... ... ... ... ...
京急富岡 -10.50 3.89 0.06 -18.09 -3.08
磯子 -10.76 2.27 0.03 -15.08 -6.22
並木北 -11.06 4.93 0.07 -21.17 -1.96
六浦 -11.13 4.24 0.07 -19.72 -3.09
並木中央 -11.40 3.92 0.06 -19.53 -4.00

133 rows × 5 columns

ちなみに、結果を表示している関数は以下。

# MCMCの結果を表示する際、パラメータ変数名ではなく、路線・駅名で表示するための関数。
def print_effectiveness(df, trace, upper_lower_station_n=None):
    # パラメータと家効果の対応表を作成する。
    beta_dict = {
        'beta__0':'切片',
        'beta__1':'駅からの距離',
        'beta__2':'築年',
        'beta__3':'部屋数',
        'beta__4':'床面積'
    }
    # パラメータと路線の対応表を作成する。
    train_dict = (df[['train','train_raw']]).drop_duplicates().set_index('train').to_dict()['train_raw']
    train_rename_dict = {}
    for key, name in train_dict.items():
        train_rename_dict.setdefault('rt__%d'%(key-1))
        train_rename_dict['rt__%d'%(key-1)] = name
    # パラメータと駅の対応表を作成する。
    station_dict = (df[['station','station_raw']]).drop_duplicates().set_index('station').to_dict()['station_raw']
    station_rename_dict = {}
    for key, name in station_dict.items():
        station_rename_dict.setdefault('rs__%d'%(key-1))
        station_rename_dict['rs__%d'%(key-1)] = name
    # traceのパラメータを家スペック・駅・路線名に置き換える
    ret_h = pm.summary(trace, varnames=['beta'])
    ret_s = pm.summary(trace, varnames=['rs'])
    ret_t = pm.summary(trace, varnames=['rt'])
    ret_h.rename(index=beta_dict, inplace=True)
    ret_t.rename(index=train_rename_dict, inplace=True)
    ret_s.rename(index=station_rename_dict, inplace=True)
    # 家スペック効果を表示する
    display(ret_h)
    # 路線効果、駅効果を表示する。
    if upper_lower_station_n is None:
        # 全部表示
        display(ret_t.sort_values('mean',ascending=False))
        display(ret_s.sort_values('mean',ascending=False))
    else:
        # 上位,下位を表示する。
        max_rows = pd.get_option('display.max_rows')
        pd.set_option('display.max_rows', upper_lower_station_n*2)
        display(ret_t.sort_values('mean',ascending=False))
        display(ret_s.sort_values('mean',ascending=False))
        pd.set_option('display.max_rows', max_rows)

5. 結論

トップでも書いたが、ADVIは圧倒的に高速だが、やはり近似であるため、MCMCの結果とは多少異なる。 ただ、傾向はつかめると思うので、検討初期の試行錯誤中はADVIを使い、ある程度モデルが固まってきたらMCMCでしっかり推定する、 というのが良いと思われる。

以上