読者です 読者をやめる 読者になる 読者になる

About connecting the dots.

statistics/machine learning adversaria.

中古マンション売買データを分析してみた(3.2) - 駅が複数の路線に所属する階層モデル

引き続き,階層モデルを実際のデータに即した形で組んでみようという試みを続けていきたいと思います.昨日の記事では,1つの駅は1つの路線にしか所属しない*1という制約を持ったモデルを作成しました.で,今日の試みは,1つの駅が複数の路線に所属できるようにモデルを拡張しましょうか,というお話です.

1つの子が複数の親に所属可能な階層モデル

さっそくですが,今回作成したのは以下のようなモデルです.1行目は昨日のものと同じですが,2行目の後半部分が総和記号を使った形になっています.これは駅の効果が固定項と,当該駅が所属する駅全てのランダム効果の総和になっているという意味です.この形をとることにより,1つの駅が複数の路線に所属することが可能になります(とりあえずここでは,路線がたくさん通っているほど,駅の魅力は線形和で増加するというシンプルなモデルを想定しています*2).

価格_{i,j,k}=b_0+b_1 距離_i+b_2 築年_i+b_3 部屋数_i+b4 床面積_i+駅_j
駅_j=r_0+\sum_{k=0}^K路線_{j,k}

モデルの記述とシミュレーション

モデル式

ということで,Stanでこれをモデリングしてみましょう.Stanコードはこちらです.前回との違いは,dataブロックの8行目がstation-train matrixとなっているように,路線とその所属駅の関係を1/0の行列で表したものになっている点です.それに合わせてparameterブロックの5行目がvector[N_T]に,modelブロックの11行目の回帰モデルの内部が行列演算の形になっています.

data {
    int<lower=1>            N;      # sample num
    int<lower=1>            M;      # independents' num
    int<lower=1>            N_T;    # train num
    int<lower=1>            N_S;    # station num
    matrix[N, M]            X;      # independents
    vector[N]               Y;      # dependent
    matrix[N_S, N_T]        ST;     # station-train matrix
    int<lower=1, upper=N_S> S[N];   # station
}
parameters {
    real             a;
    vector[M]        b;
    real             as;
    vector[N_T]      r_t;
    real             r_s[N_S];
    real<lower=0>    s;
    real<lower=0>    s_rs;
    real<lower=0>    s_rt;
}
model {
    # regresion model with random effect
    for (i in 1:N)
        Y[i] ~ normal(a+X[i]*b+r_s[S[i]], s);
    # prior distributions
    s ~ uniform(0, 1.0e+4);
    a ~ normal(39, 1.0e+4);
    for (i in 1:M)
        b[i] ~ normal(0, 1.0e+4);
    for (i in 1:N_S)
        r_s[i] ~ normal(as+ST[i]*r_t, s_rs);
    # hierarchical prior distribution
    s_rs ~ uniform(0, 1.0e+4);
    as   ~ normal(0, 1.0e+4);
    for (i in 1:N_T)
        r_t[i] ~ normal(0, s_rt);
    # 2 hierarchical prior distibution
    s_rt ~ uniform(0, 1.0e+4);
}

そしてこれをキックする側のRコードは,以下のようになります*3

# load data
d = read.delim('data/mantions.csv', header=T, sep=',')
d = na.omit(d)
attach(d)
st = read.delim('data/station_train.csv', header=F, sep=',')

# package data for stan
X  = t(rbind(distance, from, room, space))
ST = st
S  = station
Y  = price
d.stan = list(N=nrow(X),
              N_T=length(unique(train)),
              N_S=length(unique(station)),
              M=ncol(X),
              X=X,
              ST=ST,
              S=S,
              Y=Y)

5行目で駅-路線行列のstを読み込んでますが,データの中身は以下のような感じです.行が駅,列が路線で駅が路線に所属していれば1に,そうでなければ0になる行列ですね.これはPythonで前処理をして別に作っています.

http://f.st-hatena.com/images/fotolife/S/SAM/20141111/20141111220703_original.png

前処理のPythonスクリプトはこちらです.メインのモデル作成に使っているデータフレームに対して,軽く前処理して行列を作成しています.ところどころ数字がハードコーディングしてあるのは,使い捨てスクリプトだからなので,大目にみてください.

#!/usr/bin/env python
#-*-coding:utf-8-*-

import argparse
import numpy as np

# args
parser = argparse.ArgumentParser()
parser.add_argument('in_path')
parser.add_argument('out_path')


if __name__ == '__main__':
    # setup
    args = parser.parse_args()
    in_file = open(args.in_path, 'r')

    # create station-train matrix
    st = np.zeros((133, 23))
    is_header = True
    for line in in_file:
        if is_header is True:
            is_header = False
            continue
        tokens = line.split(',')
        station, train = int(tokens[11])-1, int(tokens[12])-1
        st[station][train] = 1

    # write out
    np.savetxt(args.out_path, st, fmt="%.0f", delimiter=',')

    # tear down
    in_file.close()

シミュレーションの結果

結果からいうと,これでもやはりダメでした...推定値のサマリについて,抜粋を以下に載せますが,やはり定数項と,駅の項が基本的に収束しない感じですね.まぁそれでも前回よりはだいぶマシなわけですが...

               mean se_mean      sd       2.5%        25%        50%        75%      97.5% n_eff  Rhat
a         -1549.051  68.070 137.039  -1684.644  -1636.850  -1600.390  -1539.373  -1150.424     4 2.295
b[1]         -4.784   0.001   0.098     -4.977     -4.848     -4.784     -4.718     -4.588  5001 1.000
b[2]          0.850   0.000   0.011      0.829      0.843      0.850      0.857      0.871  4491 1.000
b[3]         -0.422   0.003   0.190     -0.782     -0.554     -0.421     -0.292     -0.047  4878 1.002
b[4]          1.935   0.002   0.149      1.642      1.835      1.936      2.032      2.232  4580 1.001
as         -109.097  67.803 135.354   -510.361   -108.223    -59.499    -22.132     11.961     4 2.421
r_t[1]        2.841   0.062   4.104     -5.014      0.006      2.818      5.660     10.866  4438 1.000
r_t[2]       -0.770   0.058   3.782     -8.394     -3.210     -0.696      1.698      6.654  4208 1.000
r_t[3]        2.725   0.027   1.692     -0.631      1.602      2.733      3.893      5.989  3803 1.000
...
r_t[21]      -4.143   0.043   2.913     -9.847     -6.072     -4.144     -2.154      1.537  4646 1.000
r_t[22]      -4.474   0.037   2.474     -9.423     -6.114     -4.483     -2.814      0.415  4384 1.000
r_t[23]       2.929   0.032   2.022     -0.984      1.543      2.883      4.320      6.887  3952 1.002
r_s[1]     -111.459  67.855 135.379   -511.946   -110.823    -61.673    -24.665      9.527     4 2.424
r_s[2]     -113.088  67.878 135.432   -513.094   -112.537    -63.455    -25.973      8.144     4 2.423
r_s[3]     -103.502  67.866 135.411   -504.068   -102.333    -53.603    -16.753     17.481     4 2.423
...
r_s[131]   -107.190  67.873 135.434   -507.730   -106.498    -57.521    -20.353     13.884     4 2.423
r_s[132]   -107.515  67.853 135.364   -507.870   -106.415    -58.097    -20.576     13.855     4 2.423
r_s[133]   -111.187  67.771 135.336   -510.871   -111.277    -61.242    -24.794     10.734     4 2.419
s             7.120   0.001   0.075      6.978      7.070      7.120      7.170      7.269  4594 1.000
s_rs          4.514   0.005   0.371      3.841      4.256      4.493      4.746      5.281  4638 1.001
s_rt          6.977   0.019   1.300      4.885      6.062      6.839      7.736      9.934  4554 1.000
lp__     -12201.878   0.158   9.823 -12222.260 -12208.061 -12201.493 -12195.220 -12183.965  3888 1.000

traceplotを眺めてみても,各シミュレーションごとにめちゃくちゃな動きをしているわけで,これだと収束しないよなぁ...と思わざるを得ません.

http://f.st-hatena.com/images/fotolife/S/SAM/20141111/20141111221702_original.png
http://f.st-hatena.com/images/fotolife/S/SAM/20141111/20141111221703_original.png

これからどうするか

現状は,ここで詰まっているわけです.路線の効果を単純なランダム効果の形でモデルに突っ込むのが問題なのかなぁという気はしているので,違う形のモデルにできたらしたいのですが,さりとてどういうモデルが望ましいのかがわからず...そもそも子と親が多対多の形になる階層モデリングが間違っているんじゃないかといわれたら,それまでのような気もしますし...StanのマニュアルとかBUGS Bookとかに,こういう場合のモデル化の例とか載ってたりするんですかね...

ということで,現状暗礁に乗り上げているので,まだ今後も続くかもしれませんが,何かしらの解決策が見つからなければ次のエントリはしばらく上がらないと思います.そんなこんなでボチボチやっていきますです,はい.

*1:しかも,不動産屋が勝手に決めた(であろう)最寄駅の所属路線に従ったモデルです

*2:もちろんこれ以外にも,複数所属する路線のうちその最大値を用いる,であったり,1番目,2番目と影響力が減衰するような重み付けをして足し合わせる,といったやり方があるかなぁとは思います.

*3:コードの全体は例のごとくgithubにあげてあるので,こちらを参照してください.