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

About connecting the dots.

statistics/machine learning adversaria.

中古マンション売買データを分析してみた(2.1) - MCMCで線形モデルを推定するためのパラメタ設定はどうすればいいか

今回のエントリはテクニカルな問題のお話のみで,住宅価格についての蘊蓄はいっさい含まれておりません.統計モデリングに興味のある方のみお読みください.

問題

前回のモデルでは,単純に群間の差を駅,路線の2層の群を使ってモデリングしました.考察のところでも書いたように,実のところあのモデルは余りよろしくないわけです.なにかといったら,前々回モデリングで使っており,有用であることがわかっている駅からの距離,築年数という変数を使っておらず,そのため推定結果にバイアスがかかっていました*1.じゃあやれよ,という話になりますが,実際には,やってみた上でパラメタ収束が劣悪だったため,モデルに組み込むのをあきらめたという経緯があります.

しかし,普通にlm()で導出した線形回帰でも,十分に有意な値が出ているし,分布形状をみてもそこまで当てはまりが悪いとは考えにくいわけです(もちろん厳密には正規分布近似が妥当か? というような形状であることは否定しませんが...).

http://f.st-hatena.com/images/fotolife/S/SAM/20140531/20140531110237_original.jpg

分析結果

lm()による最小二乗法推定

lm()を使って推定した結果が以下になります.全体的にモデルの当てはまりはそれなりで,まぁ悪くない数字なのではないかと思います,

> summary(lm(price~space+room+distance+from, data=d))

Call:
lm(formula = price ~ space + room + distance + from, data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-30.449  -6.879  -0.793   5.841  75.247 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.873e+03  2.667e+01  -70.23   <2e-16 ***
space        2.796e+00  1.987e-01   14.07   <2e-16 ***
room        -2.608e+00  2.510e-01  -10.39   <2e-16 ***
distance    -4.230e+00  1.200e-01  -35.24   <2e-16 ***
from         9.623e-01  1.343e-02   71.63   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.06 on 4821 degrees of freedom
Multiple R-squared:  0.6451,	Adjusted R-squared:  0.6448 
F-statistic:  2190 on 4 and 4821 DF,  p-value: < 2.2e-16

JAGSによるMCMC推定

1500 iteration, 100 burn-in, 3 thin*2

さて,これを同様にJAGSを使って推定してみます.もちろんこんな単純なモデルをMCMC推定する意味は全くないのですが,この先に一般化線形混合モデル,階層ベイズモデルに展開することを考えているため,まずはシンプルなモデルをJAGSで組むことにしたわけです.コードは以下の通りです*3

# mcmc modeling
## model description
mantion_model <- function() {
  # main model
  for (i in 1:N) {
    PRICE[i] ~ dnorm(b+bD*DISTANCE[i]+bF*FROM[i]+bR*ROOM[i]+bS*SPACE[i], tau[1])
  }
  # prior distribution
  b  ~ dnorm(0.0, 1.0e-6);
  bD ~ dnorm(0.0, tau[2]);
  bF ~ dnorm(0.0, tau[3]);
  bR ~ dnorm(0.0, tau[3]);
  bS ~ dnorm(0.0, tau[4]);
  for (i in 1:N.tau) {
    tau[i] <- 1 / (sigma[i] * sigma[i])
    sigma[i] ~ dunif(0.0, 1.0e+6)
  }
}

# モデルの書き出し先
file.bugs <- file.path('script/mantion.bug')
R2WinBUGS::write.model(mantion_model, file.bugs)

# data for JAGS
list.data <- list(
  PRICE      = d$price,
  DISTANCE   = d$distance,
  FROM       = d$from,
  ROOM       = d$room,
  SPACE      = d$space,
  N          = nrow(d),
  N.tau      = 4
)

# initial value for JAGS
inits <- list(
  b     = rnorm(1, 39.22, 16.88),
  bD    = rnorm(1, 0, 100),  
  bF    = rnorm(1, 0, 100),
  bR    = rnorm(1, 0, 100),
  bS    = rnorm(1, 0, 100),
  sigma = runif(4, 0, 100)
)

# estimate parameter
params <- c('b', 'bD', 'bF', 'bR', 'bS', 'sigma')

# model definition
model <- jags.model(
  file    = file.bugs,
  data    = list.data,
  inits   = list(inits, inits, inits),
  n.chain = 3
)

# set MCMC burn-in cycles
update(model, 100)

# conduct MCMC silumation results
post.list <- coda.samples(
  model,
  params,
  thin = 3,
  n.iter = 1500
)
plot(post.list)

推定結果のプロットを行うと,全く推定結果が収束していないのが見て取れるかと思います.特に切片と築年が酷いことになっています.ちゃんと理解できていないのですが,これってMCMCイテレーションをもっとずっと増やせば安定するんですかね? 100万回とか...

http://f.st-hatena.com/images/fotolife/S/SAM/20140531/20140531114551_original.jpg
http://f.st-hatena.com/images/fotolife/S/SAM/20140531/20140531114552_original.jpg
http://f.st-hatena.com/images/fotolife/S/SAM/20140531/20140531114553_original.jpg

実際に各パラメタ推定におけるR-hatを確認して,収束診断をしてみましょう.JAGSではR-hatを確認できないのですが,@beroberoさんに教えていただいたやり方で,WinBUGSパッケージのオブジェクトに変換してR-hatを確認することができます.やはりbとbFの収束が1.5と悪いことが見て取れます.

推定された各回帰係数も,最小二乗法とは大きく異なっているので,いろいろとアレな感じがしますね.とりあえず,もう少し収束っぽくなるまで計算回し続けてみますか...

# convert to BUGS object
mcmc.list2bugs <- function(mcmc.list) {
  b1 <- mcmc.list[[1]]
  m1 <- as.matrix(b1)
  mall <- matrix(numeric(0), 0, ncol(m1))
  n.chains <- length(mcmc.list)
  for (i in 1:n.chains) {
    mall <- rbind(mall, as.matrix(mcmc.list[[i]]))
  }
  sims.array <- array(mall, dim = c(nrow(m1), n.chains, ncol(m1)))
  dimnames(sims.array) <- list(NULL, NULL, colnames(m1))
  mcpar <- attr(b1, "mcpar")
  as.bugs.array(
    sims.array = sims.array,
    model.file = NULL,
    program = NULL,
    DIC = FALSE,
    DICOutput = NULL,
    n.iter = mcpar[2],
    n.burnin = mcpar[1] - mcpar[3],
    n.thin = mcpar[3]
  )
}
> mcmc.list2bugs(post.list)

 3 chains, each with 2600 iterations (first 1100 discarded), n.thin = 3
 n.sims = 1500 iterations saved
           mean      sd   2.5%    25%    50%    75%   97.5% Rhat n.eff
b        -211.5    26.8 -264.8 -232.6 -210.0 -191.1  -161.7  1.5     7
bD         -5.5     0.2   -5.8   -5.6   -5.5   -5.4    -5.1  1.0   110
bF          0.1     0.0    0.1    0.1    0.1    0.1     0.2  1.5     7
bR         -4.2     0.3   -4.8   -4.4   -4.2   -4.0    -3.5  1.0    76
bS          6.6     0.3    6.0    6.4    6.6    6.8     7.1  1.0   110
sigma[1]   13.5     0.2   13.2   13.4   13.5   13.6    13.8  1.1    26
sigma[2] 2568.7  4074.1    5.0   49.4  569.9 2995.7 15165.5  1.0   250
sigma[3]   14.9    27.4    1.8    3.6    6.1   14.0    95.2  1.0   460
sigma[4] 5500.0 10019.8    6.4   98.0 1079.4 6296.9 41095.6  1.1    50

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
100000 iteration, 100 burn-in, 3 thin

先ほどのシミュレーションだと,bとbFの傾きが変化し続けている状態で,明らかにまだ収束に至っていなさそうだったので,イテレーション回数を増やして再計算することにしました.3回まわしで各10万イテレーションにして計算しなおしたところ,今度はおおむね収束したといえそうな形になりました.値もおおむねlm()のときと一致しており,妥当な形になったといえそうです.実際に収束に至るまでには,だいたい2万イテレーションぐらいあれば良かったみたいですね.まぁイテレーション回し終わるのに6時間くらいかかりましたが...

> mcmc.list2bugs(post.list)

 3 chains, each with 1001102 iterations (first 1100 discarded), n.thin = 3
 n.sims = 1000002 iterations saved
            mean      sd    2.5%     25%     50%     75%    97.5% Rhat n.eff
b        -1830.4   159.7 -1915.0 -1884.2 -1867.2 -1846.9  -1267.3    1 11000
bD          -4.3     0.2    -4.7    -4.3    -4.2    -4.2     -4.0    1 20000
bF           0.9     0.1     0.7     0.9     1.0     1.0      1.0    1 22000
bR          -2.6     0.3    -3.3    -2.8    -2.6    -2.4     -2.1    1 94000
bS           2.9     0.4     2.4     2.7     2.8     3.0      4.2    1 10000
sigma[1]    10.1     0.2     9.9    10.0    10.1    10.1     10.6    1 60000
sigma[2] 26003.6 56848.5     4.2    62.6   970.6 15224.3 216815.1    1   400
sigma[3]    13.2    47.1     1.2     2.4     4.1     8.7     79.7    1  8300
sigma[4] 24431.4 46458.6     3.0    55.7  1129.1 23171.6 161349.9    1    86

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

http://f.st-hatena.com/images/fotolife/S/SAM/20140531/20140531193151_original.jpg
http://f.st-hatena.com/images/fotolife/S/SAM/20140531/20140531193150_original.jpg
http://f.st-hatena.com/images/fotolife/S/SAM/20140531/20140531193152_original.jpg

80000 iteration, 20000 bunr-in, 40 thin

ということで,今度は値が変化している状態の不安定なパラメタをburn-inに入れてしまって,さらに間引きも40にして計算を軽くしてみたいと思います.そうすると,上の場合とそれほどは変わらない値が得られ,かつ計算時間も2時間程度とだいぶ高速化されました.ただbのR-hatがまだ1.1というところをみて,もう少し間引き数を減らすかステップ数を増やすかの対応が必要そうです.

> mcmc.list2bugs(post.list)

 3 chains, each with 101000 iterations (first 21000 discarded), n.thin = 40
 n.sims = 6000 iterations saved
             mean       sd    2.5%      25%      50%      75%    97.5% Rhat n.eff
b         -1712.2     91.0 -1829.2  -1790.2  -1739.4  -1653.5  -1517.0  1.1    50
bD           -4.4      0.1    -4.6     -4.4     -4.4     -4.3     -4.1  1.0   220
bF            0.9      0.0     0.8      0.9      0.9      0.9      0.9  1.0    52
bR           -2.8      0.3    -3.3     -2.9     -2.8     -2.6     -2.3  1.0   340
bS            3.2      0.3     2.7      3.0      3.1      3.3      3.8  1.0    81
sigma[1]     10.1      0.1     9.9     10.0     10.1     10.2     10.3  1.0   440
sigma[2] 493966.3 289142.0 20444.8 244921.5 488970.6 745387.6 973813.5  1.0  6000
sigma[3] 496396.3 284594.6 25977.7 252167.8 493177.6 739757.1 974007.0  1.0  6000
sigma[4] 498123.5 288435.7 22577.8 245519.2 502556.0 749505.3 972328.5  1.0  5100
sigma[5] 505126.8 289866.8 25836.7 251756.4 507867.7 760266.6 976506.2  1.0  6000

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

ということで

自分の理解があやふやなのに問題があるのはわかっているわけですが,とりあえず試行錯誤の過程を公開しておくことにします.引き続き精進していきます.最後に,アホとしてこちらのメッセージを胸に刻んでおきたいと思います.高々線形回帰するだけで6時間計算にかけてるようじゃね...Stan使わないとこれ以上の処理は無理なのでは感.

*1:より具体的には,みなとみらいと武蔵小杉のあの突出っぷりは,は,両者が新興住宅地かつ駅近タワーマンション乱立地帯ということが大きく影響していると思われ,本来は駅からの距離と築年数をモデルに入れこんで,その影響をのぞいた上での,両駅のブランド力を測定したかったのです.

*2:iterationは,MCMCシミュレーションを実施するステップ数,burn-inはパラメタが安定するまでの初期ステップを無視するステップ数,thinは計算結果をハンドリングしやすくするための間引き間隔を表します.

*3:例のごとく全体はgithubを更新済みです.