Data Science by R and Python

統計学を、広く、深く、わかりやすく。

pythonでstepwise regression用の関数を作っておく

はじめに

ほんと、久々の更新になってしまいました。。。
いまだに月間で1000PVほど見られているようでとてもありがたく思いますm(_ _)m

最近も変わらず因果推論の研究を中心に行っておりますが、それ関連の内容はまた機会をみてblogで書いていければと思っています。
また先日、twitterで公開したこちらのスライドもたくさんの方に見ていただけまして、コメントも頂けたりして、とても嬉しく、励みになっています。
speakerdeck.com

また、少しずつではありますが更新いたしますので、たまに覗いていただければ嬉しいです。
では、本題にまいります。

今回の更新

とはいっても、今日の更新は、大した内容ではなく、pythonでstepwise regressionの関数で自分がほしいものがないので、つくりましたという内容です。

Stepwise Regressionについて

特に、回帰モデルでモデルに含める変数をある基準に基づいて選択するという方法です。
有名なのは、AICを用いるものですが、今回はBIC(Bayesian Information Criterion)で実装します。

pythonのstepwise regressionについて

一応、stepwise regression用のモジュールがあるのですが、p-valueをみて変数を選択していて、ちょっと「うーん?」という感じ。
あと、調べると他にもあるのですが、statmodelsに依存した形のものが多かったので、できればそのあたりも依存しないで動かしたくて、結局自作に至りました😣
このあたりは、Rはほんとしっかりしていて、step functionなんかを使うとサクッとできるあたりは嬉しいところですね。

本題 - 関数 -

とりあえずAICBICを基準に選択できればOKということと、この関数を別の関数の中で使用するために、返り値は回帰係数であってほしいという個人的な理由が存分に入っております!
Rっぽい出力はまた時間があれば作成しますが、やる気が風前の灯なのでこれくらいで妥協いたします...ご勘弁ください!!

# numpy
import numpy as np

#Akaike Information Criterion
def AIC_gaussian(x,y):
    y_pred = np.dot(x,np.linalg.solve(np.dot(x.T,x),np.dot(x.T,y)))
    resid = y - y_pred #residuals
    rss = np.dot(resid.T,resid) #residual sum of square
    AIC = len(y)*np.log(rss/len(y)) + (x.shape[1])*2
    return AIC

#Bayesian Information Criterion
def BIC_gaussian(x,y):
    y_pred = np.dot(x,np.linalg.solve(np.dot(x.T,x),np.dot(x.T,y)))
    resid = y - y_pred #residuals
    rss = np.dot(resid.T,resid) #residual sum of square
    BIC = len(y)*np.log(rss/len(y)) + (x.shape[1])*np.log(len(y))
    return BIC

def forward_stepwise_result(x,y):
    beta = np.linalg.solve(np.dot(x.T,x),np.dot(x.T,y))
    return beta

# stepwise selection
def forward_stepwise(x,y,method='BIC',verbose=True):
    x_tmp = np.concatenate([np.reshape(np.ones(1000),(1000,1)),x],axis=1) #xにinterceptを追加
    initial_list = [0]
    included = list(initial_list)
    while True:
        changed=False
        excluded = list(set(np.arange(x_tmp.shape[1]))-set(included))
        BIC_res = np.zeros(len(excluded))
        base_BIC = BIC_gaussian(x_tmp[:,included],y)
        print('Baseline model with BIC: {:}'.format(base_BIC))
        j = 0
        for new_column in excluded:
            BIC_res[j] = BIC_gaussian(x_tmp[:,included+[new_column]],y)
            print('Add :{:30} with BIC: {:}'.format(excluded[j],BIC_res[j]))
            j = j+1
        if BIC_res.min() < base_BIC:
            best_feature = excluded[BIC_res.argmin()]
            included.append(best_feature)
            changed=True
            if verbose:
                print('Finally Add  {:30} with BIC {:}'.format(best_feature, BIC_res.min()))

        if not changed:
            if verbose:
                print('Any variables does not added, stop forward stepwise regression')
            break
            #end while
    
    #final result の回帰係数を返す
    beta = np.reshape(np.zeros(x_tmp.shape[1]),(x_tmp.shape[1],1))
    beta[included] = forward_stepwise_result(x_tmp[:,included],y)
    
    return beta

#結果
forward_stepwise(x,y,method='BIC',verbose=True)

#Baseline model with BIC: [[7214.43429935]]
#Add :                             1 with BIC: 6564.464745289761
#Add :                             2 with BIC: 7062.695910602746
#Add :                             3 with BIC: 7220.4285978234575
#Add :                             4 with BIC: 6830.525994981331
#Finally Add                               1 with BIC 6564.464745289761
#Baseline model with BIC: [[6564.46474529]]
#Add :                             2 with BIC: 6260.368109310093
#Add :                             3 with BIC: 6571.08187071489
#Add :                             4 with BIC: 5708.315703761152
#Finally Add                               4 with BIC 5708.315703761152
#Baseline model with BIC: [[5708.31570376]]
#Add :                             2 with BIC: 5713.846539819814
#Add :                             3 with BIC: 5710.261885909516
#Any variables does not added, stop forward stepwise regression
#array([[21.86610567],
#       [42.05186786],
#       [ 0.        ],
#       [ 0.        ],
#       [ 0.35454755]])

参考にしたものなど

ベイズ情報量規準 - Wikipedia

  • p-valueを使ってるので、欲しいものではなかった....

github.com

  • 多変量モデルの選択

www.asakura.co.jp


それでは、ハッピーメリークリスマス!!

因果効果の推定!Rで実践 - 傾向スコア,マッチング,IPW推定量 -

はじめに

先日,岩波データサイエンスvol.3が発売されました.私も,「傾向スコアを用いたバント効果の推定−−ノーアウト1塁のバントは,得点確率を有意に高めるか!?」と題した記事を寄稿させていただきました.興味がある方は,是非読んでください&感想・コメントなどもTwitterなどで教えてください.

岩波データサイエンス Vol.3

岩波データサイエンス Vol.3

今回のテーマ

今回は,上記の記事に関連して,上記の本で用いているIPW推定量の推定方法を実際に実践してみます.
説明していく過程で,読んでくださる方に,
・2つの群の平均の比較で因果効果をきちんと推定できない例
・傾向スコアの推定方法(ロジスティック回帰を今回は用います)とその周辺知識
・マッチングの計算方法
・IPW推定量の計算方法
をご紹介したいと思います.ただ,詳しく書きすぎると長くなってしまうので,詳細な部分については参考文献をあげるにとどめ,本の内容にかぶってしまうところは省きます(買って読んでください!お願いします).では,始めましょう.

データセットをつくる

今回は,適当な処置d=0,1が結果変数Yに与える影響を推定するという問題を考えます(d=1のときは処置を行い,d=0のときは処置を行わなかったとします).ただ,普通の例では面白くないので,
・実際には「処置dはYを減少させる効果がある」のに「単純な2群の比較をするとYを増加させる影響がある」という結論を導いてしまうという例を作成します.

こういう例は,実は比較的簡単に生成できます.まずは,共変量を作りましょう.ここでは,適当で構いません.とりあえず正規分布に従う共変量を4つほど作っておきます(先に必要な関数を読み込んでおきます.IPW計算用の関数と,ロジット変換などを手軽にするための関数です.)

> ######functions######
> expit = function(x){
+ 	rr = exp(x)/(1+exp(x))
+ 	return(rr)
+ }
> 
> logit = function(x){
+ 	rr = log(x/(1-x))
+ 	return(rr)
+ }
> 
> IPW = function(n,p,d,y){
+ 	numerator = sum(d*y/p)
+ 	denominator = sum(d/p)
+ 	r = numerator/denominator
+ 	return(r)
+ }
>#############
> set.seed(20160619)
> 
> #------------------------------
> # 標本数
> #------------------------------
> 
> n = 800
> 
> #------------------------------
> # 共変量
> #------------------------------
> 
> x1 = rnorm(n)
> x2 = rnorm(n)
> x3 = rnorm(n)
> x4 = rnorm(n)

ここからがポイントです.この共変量を使って,傾向スコアの真値に対するモデルを作成します.改めて述べておきますが,傾向スコアというのは,共変量が観測されたもとで,処置を受ける確率 p(d=1|x)です.今回は,この確率に対して,以下のモデルを仮定します.

e =  p(d=1|x) = \cfrac{\exp(-0.8\times x_1 + 0.8 \times x_2 + 0.25 \times x_3 - 0.2 \times x_4)}{1+\exp(-0.8\times x_1 + 0.8 \times x_2 + 0.25 \times x_3 - 0.2 \times x_4)}

ここで,この式を見てもらうと, x_{2}, x_{3}が大きいほど処置を受ける確率が高く, x_{1}, x_{4}が小さいほど処置を受ける確率が高くなるということがわかります(と,作っただけですが).そして,この傾向スコアの真値に基づいて割り付け変数 dを,ベルヌーイ分布から発生させます.

> #------------------------------
> # 傾向スコアの真値
> #------------------------------
> 
> p = expit(-0.8*x1 + 0.8*x2 + 0.25*x3 - 0.2*x4)
> 
> #------------------------------
> # 各標本の割り付け変数(d=1,0)
> #------------------------------
> 
> d = rbinom(n,1,p)

次に,共変量と割り付け変数を用いて結果変数yを作ります.ここで,傾向スコアが大きいときほど,Yが大きな値をとるようにします.具体的には, x_{2}, x_{3}が大きいほど処置を受ける確率が高く, x_{1}, x_{4}が小さいほど処置を受ける確率が高くなるわけですから,yをつくるときには, x_{2}, x_{3}の係数を正に, x_{1}, x_{4}の係数を負にするとよいです.こうすると,傾向スコアが高いものほど,yが大きくなりますから,処置の効果がマイナスでも,係数さえ調整すれば,処置群と対照群を平均で比較すると,処置群のほうがyの平均は大きくなります.

ここでは,yに対するモデルとして,
 y = 80 + (-7) \times d - 7 \times x_{1} + 6\times x_{2} + 6 \times x_3 - 4 \times x_{4} + \mbox{誤差}
とします.

> #------------------------------
> # 結果変数 d=1なら結果は7下がる
> #------------------------------
> 
> y = 80 + (-7)*d - 7*x1 + 6*x2 + 6*x3 - 4*x4 + rnorm(n)
> 


最後に,観測データを作りますが,実際に観測できる変数は x_{1}, x_{2},x_{3}, x_{4}と,d, yですから(傾向スコアの真値は,観測されません),次のような観測データ行列を作っておきます.

> #------------------------------
> # 観測データ
> #------------------------------
> 
> data = cbind(x1,x2,x3,x4,d,y) 
> head(data)
             x1         x2         x3         x4 d        y
[1,] -1.2163990 -1.9305651 -0.6880511  0.1844445 0 72.59452
[2,]  0.2341120  0.1233155  0.9854147 -1.6495960 1 84.68886
[3,]  0.4495428  0.2482259  0.3206955 -0.2648300 1 75.93403
[4,]  0.6136386 -0.1667101 -0.9992488 -0.2649914 0 70.71651
[5,] -0.9403283  0.2755147  1.4383773  0.2877559 1 89.74142
[6,] -0.7335436  0.3948032 -1.3212992  0.6402797 1 70.53397

単純な2群の比較

では,実際に観測されたデータでd=1の群(処置群)とd=0の群(対照群)の結果変数を比較します.

> #------------------------------
> # 群ごとの平均の差
> #------------------------------
> 
> mean(y[d==1]) - mean(y[d==0])
[1] 3.256249

結果は正になりました.つまり,2群の平均を比較するということを行うと,処置の効果によってYの値は上昇するという結論を導くことになります.しかし,先ほどのyに対する式を確認しますと,
 y = 80 + (-7) \times d - 7 \times x_{1} + 6\times x_{2} + 6 \times x_3 - 4 \times x_{4} + \mbox{誤差}
dの係数はマイナスなので,d=1である,つまり処置を受けた場合には,Yの値は減少するはずなのです.つまり,2群の比較によって得られる結果は,処置の因果効果を正しく推定できていません.

そこで,傾向スコアの登場です.ここからは,傾向スコアの推定と,それを用いた因果効果の推定法について述べます.

傾向スコアの推定と,傾向スコアを利用した推定量

まず,傾向スコア推定に焦点を当てます.傾向スコアとは,共変量が与えられた元で,処置群に割り付けられる確率でした.少し,わかりにくいかもしれませんが,言い換えると「共変量xの標本は,確率e(x)で,表が出るコインを振って,表なら処置群に,裏なら対照群に入る」ということです.このような問題は分布の言葉で書くと,ベルヌーイ分布を用いて,次のように書けます.
 (d|x) \sim \mbox{bernoulli}(e(x_{i}))

気づいた方もいるかもしれませんが,これは一般化線形モデルの形に似ています.しかし,e(x)が共変量xのどのような関数型で表現できるかは実際にはわからないので,このままでは,以下のようになります
e(x) = f(x)

ただ,関数型fがわからないと,傾向スコアを推定できません,そこで,ここでは,ロジスティック回帰モデルを選択するということにします.つまり,以下のように関数の形に仮定を置くことにします.
 \log\cfrac{e(x)}{1-e(x)} = x'\beta

荒い言い方になりますが,一般的には傾向スコアが,どんな形でxに依存するかはわからないので,試す必要がある場合もあります.また,どの変数xに依存しているのかもわからないという問題もあります.その場合も,とりあえず関係のありそうな変数を入れて,また対数変換,交互作用,指数変換などを用いて,傾向スコアのもっともらしいモデルを探します.

その際には,モデルのただしさを何で判断するの?という問題があります.これについては各分野の人が,「目安」を用意していますが,この基準さえ満たしておけば良いという判断基準はないので,あくまで目安であることに注意した方がよいと思います(実際,うまく傾向スコアを予測できる場合もありますが,そうではない場合もあります.データや分野に依存するのではないでしょうか?).

では,観測されたデータを用いて,ロジスティック回帰を当てはめて傾向スコアを推定します.

> #------------------------------
> # 傾向スコアの推定
> #------------------------------
> 
> #観測データから傾向スコアの推定を行う.
> #傾向スコアは共変量が与えられたもとでの,d=1になる確率P(d=1|x).
> #ただ,何もモデルを仮定しないままでは推定できないので,ここではロジスティック回帰を使う
> 
> #割り付け結果を,共変量で回帰する
> res = glm(d~x1+x2+x3+x4,data=as.data.frame(data),family=binomial(link="logit"))
> 
> summary(res)
Call:
glm(formula = d ~ x1 + x2 + x3 + x4, family = binomial(link = "logit"), 
    data = as.data.frame(data))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.3775  -0.9835   0.4383   0.9341   2.5166  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.15718    0.08026   1.958   0.0502 .  
x1          -0.71110    0.08596  -8.272  < 2e-16 ***
x2           0.81237    0.09085   8.941  < 2e-16 ***
x3           0.36507    0.08447   4.322 1.55e-05 ***
x4          -0.33182    0.08230  -4.032 5.53e-05 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1106.39  on 799  degrees of freedom
Residual deviance:  916.61  on 795  degrees of freedom
AIC: 926.61

Number of Fisher Scoring iterations: 4

> #ここでの傾向スコアは,ロジスティック回帰を行った時の確率の推定値(当てはめ値)
> hat.e = res$fitted
> head(hat.e)
        1         2         3         4         5         6 
0.6784036 0.8307323 0.7208275 0.7212536 0.6900221 0.4944932 
[f:id:tomoshige_n:20160619031457p:plain]

これで,傾向スコアはhat.eとして推定されたことになります.ここで,傾向スコアの推定値は次の通りです.\hat\betaはロジスティック回帰で推定された各共変量の回帰係数です.
 \hat{e}_{i}(x_{i}) = \cfrac{\exp(x_{i}'\hat{\beta})}{1+\exp(x_{i}'\hat{\beta})} .

また,最初に設定した真のモデルがロジスティック回帰モデルの平均構造でしたから,そこから発生したデータに同じモデルを当てはめていますので,推定された回帰係数は,もちろん最初に設定した真のモデルの回帰係数に近くなります.

ここからは,処置の平均処置効果(処置の因果効果)を求めていきましょう.

マッチングによる推定量の構成

マッチングは,傾向スコアが同じ標本の,処置群の標本平均と,対照群の標本平均を計算して,その差を計算して,計算結果を群ごとに傾向スコアについて平均することによって,因果効果を推定するというのが基本の考え方です.では,なぜ傾向スコアが同じ標本同士を比較するのでしょうか?マッチングの本質はここにあります.

この理由を理解するためには,そもそもの因果効果を解析する際の仮定を思い出す必要があります.それがなんだったかというと「強く無視できる割り付け」というものでした.これは,「割り付けに関係する共変量をすべて固定してしまえば,割り付け変数dと,潜在結果変数の組(y_{1},y_{0}が独立になるというものです.この結果は,極めて重要な意味を持ちます.次の式をみてください.

 E(y_{1}-y_{0}) = E\left( E(y_{1} - y_{0} | x) \right) = E\left( E(y_{1} | x)\right) - E\left( E(y_{0} | x) \right) = E\left( E(y_{1} | d =1, x)\right) - E\left( E(y_{0} | x, d=0) \right)

この式の左辺は,因果効果です.最後の等号は強く無視できる割り付けの仮定から成立します.この式の意味は,処置群と対照群で同じ共変量xを持つ標本の平均を計算した後に,その差を計算するという操作を,すべてのxで平均することができれば,平均処置効果は求めることができるという意味です.もっとわかりやすく(厳密さを欠いて言えば),「各共変量の値について,その値を持つ標本を処置群・対照群から選び出し,その標本の結果変数の差を計算するという操作を,すべての共変量の取りうる値について平均化すればいいんだ」と言っています.

これは,つまり,どちらの群に割り付けるのかを決める共変量xの値が等しい標本ごとの差を計算せよ!と言っています.つまり,その共変量を固定した箇所で見れば,標本はランダムに割り付けられていることになります.なので,これが処置の効果になるんです.しかし,すべての共変量で揃えるのは大変です...が,ここからがとても大事な結果です.実は上の式は,強く無視できる割り付けが成立する元では,xのところをe(x)(傾向スコア)に変えても成立するということが示されています(Rosembaum, and Rubin 1983).すなわち,

 E(y_{1}-y_{0}) = E\left( E(y_{1} | e(x), d =1)\right) - E\left( E(y_{0} | e(x), d=0) \right)

これはすなわち,傾向スコアが等しいところでの群平均の差を,傾向スコアで足しあわせて平均かしたものが処置の効果になると言っているのと同じことですね!そして,この考え方を利用するのがマッチングなわけです.

ただし,傾向スコアの値がぴったり同じという標本を見つけるのは難しいわけですね.そこで,傾向スコアを適当な離散区間に区切ります.今回は0.05幅に区切った結果を示します.すると,結果は-7.221となり,見事に因果効果-7に近い値が得られましたね!

しかし,ここで区間幅の決め方をどうするのかが恣意的になっていることが,ちょっと厄介ですね.実際,区切りを荒くしてみると結果は変わりますし,細かくしすぎると群ごとの平均が取れなくなるので,結果は不安定になります.これは困ったということで,別の方法が次のIPW推定量です.
(*興味がある方は,他の区切り方も試してください.)

>#----------------------------------
># 傾向スコアマッチングによる因果効果推定
>#----------------------------------
> interval = seq(0,1,0.05) #machingを行う時の傾向スコアの幅を指定する
> maching_res = rep(0,10)
> 
> for(i in 1:(length(interval)-1)){
+ 	temp0 = y[(d==0)&((interval[i]<hat.e)&(hat.e<interval[i+1]))]
+ 	temp1 = y[(d==1)&((interval[i]<hat.e)&(hat.e<interval[i+1]))]
+ 	maching_res[i] = mean(temp1) - mean(temp0)
+ }
> 
> maching_res
 [1] -5.1711970        NaN -0.3080318 -4.2663056 -8.1488541 -7.9177035 -5.4518013
 [8] -6.9564249 -6.8824864 -7.6636229 -7.1455726 -7.2535686 -7.4277173 -5.7344403
[15] -6.7976385 -8.0770380 -6.7897113 -6.7753327 -8.8878558        NaN
> 
> mean(maching_res,na.rm=T)
[1] -6.536406

IPW推定量へと

最後にIPW推定量についてです.IPW推定量は,その定義から明らかなように,傾向スコアの逆数で重みをつけて因果効果を推定します.これにはどんな意味があるのでしょうか?また,1つの標本を考えます.この標本の共変量はxであるとします.このとき,この標本が処置群に行く確率は80%で,対照群に行くか確率は20%であるとします.さらに,ここでは簡単のため処置群と対照群の人数はそれぞれ1000人であるとします.つまり,100人の同じ共変量xの人がいれば,期待値の意味で,80人は処置群へ,20人は対照群へ行くわけです.そうすると,共変量すべての標本について和をとるとき,処置群では,共変量xを持つ人の結果変数は80/1000 = 8%考慮され,対照群では20/1000=2%考慮されます.

何が言いたいかというと,共変量xを持つ人の1000人中に占める割合が異なってしまうために,両軍でそれぞれ平均を計算すると共変量のせいで差が生まれます.実際,結果変数がxに依存している場合には,当たり前ですが,差が生まれます.

今回の共変量に依存する割り付けというのは次の問題のように捉えてもいいと思います.

今80人の子供がいます.身長は150cmか,140cmの人しかおらず,150cmの人の体重は50kg, 140cmの人の体重は40kgです.そして,50kgの人は確率0.8でクラス1へ,確率0.2でクラス2に入ります.一方で,40kgの人は確率0.8でクラス2へ,確率0.2でクラス1へ入ります.そして,体重が50kgと40kgの人はそれぞれ40人です.このとき,クラス1とクラス2の平均身長を比較すると異なる理由はなぜでしょうか?もしランダムに割り付けられたとすると,クラスの平均身長はどうなるでしょうか?

こう考えると,これは体重という共変量が,割り付け(クラス分け)と,結果変数(身長)に影響をするから,クラス1とクラス2の平均身長が異なるということになります.では,補正するためにはというと,クラス1では40kgの人に1/0.2の重みで身長を評価し,50kgの人を1/0.8で評価します.これは,割り付け確率の逆確率(つまり傾向スコアの逆重み付けです)です.一方で,クラス2では50kgの人に1/0.2の重みで身長を評価し,40kgの人を1/0.8で評価します.すると,クラスの平均身長の差の平均は0になることがわかります.なぜなら,最初の割り付けにおいて,身長150cmの人は0.8の確率でクラス1に入りますから,平均すると40*0.8 = 32人がクラス1に入ることになります.そう考えると,クラス1とクラス2の調整済み平均は
クラス1:(150/0.8×32 + 140/0.2×8)/40 = 145
クラス2:(150/0.2×8 + 140/0.8×32)/40 = 145
となり,割り付けをきちんと戻せば,平均に差がないことがわかります.

これと同じことを考えているのが,傾向スコアによるIPW推定量ということになります.つまり,傾向スコアの逆数による重み付けというのは,この共変量によって影響される人数比の偏りを補正する役目(つまりは,ランダムな割り付けにするための処置!!)があるわけですね.さて,この考え方からIPW推定量は,ここで書かれているようになります.

smrmkt.hatenablog.jp

> #-----------------------------
> # IPW推定量の計算
> #-----------------------------
> 
> #傾向スコアが求まったらIPW推定量を計算します.
> IPW.Treated = IPW(n,hat.e,d,y)
> IPW.Untreated = IPW(n,1-hat.e,1-d,y)
> 
> #因果効果
> IPW.Treated - IPW.Untreated
[1] -6.953027

近い値が得られました!めでたし!(ただ,set.seedを変えると...ちょっとざんない結果になることもありますので,ご了承を!ただ,それは推定量のばらつきなので,問題ありません.)

おわりに

ここまで,長い道のりをお付き合いくださりありがとうございました.
ここまでなんと,11,000字もあり,簡単な論文程度の文字数になってしまいました.
8時間もかかってしまいました.外明るくなってきてる気がする...
傾向スコアがどういうもので,IPW推定量・因果効果とは何かについて少しでも助けになれば幸いです.

岩波データサイエンスvol.3を片手に持って,
一つ一つのRコードの意味を岩波データサイエンスで確かめながら,
読んでいただけれいれば,さらに嬉しく思う次第です.

野球データの解析について

岩波データサイエンスに寄稿させていただいた野球データの解析では,
上記の議論に加えて,平滑化スプラインなど発展的な手法を用いて解析しています.
この辺りについてはサポートページに書かせていただきましたので参照ください.
ただし,傾向スコアや,IPW推定量の計算はここに書かせていただいた方法と大きくは変わりません.

・サポートページ
【特集】因果推論 ― 実世界のデータから因果を読む - 岩波データサイエンス

そして最後に,宣伝.←しつこくてごめんなさい...

まずは,岩波データサイエンスvol.3!因果推論の素晴らしい先生方がわかりやすく,因果推論にまつわる話を網羅してくださっています.是非手にとって,内容をごらんください!.因果推論のもやもやっとしたところが解消されるかと思います!

岩波データサイエンス Vol.3

岩波データサイエンス Vol.3

それから,因果推論について詳しく知りたい方,式の展開が気になる方は,星野先生の「調査観察データの統計科学」を手に取ることをお勧めします.基礎から,発展まで数学的な記述を含めて網羅的に書かれております.

推定量の一致性(大数の法則)の感覚的に理解をする

久々の更新です.

今日,シミュレーションで久々に推定量の一致性(大数の法則)を可視化してみたので,せっかくなのでブログに載せておこうと思います.

データを発生する確率変数 {X}の分散を {\sigma^{2}}とします.このとき,確率変数{X}から発生させられたサイズの {N}標本から,真の分散 {\sigma^{2}}を推定することを考えます.このとき, {\sigma^{2}}に対する不偏な推定量(期待値の意味で, {\sigma^{2}}に一致する推定量)は,


\begin{align}
 \displaystyle \hat\sigma^{2} = \frac{1}{N-1}\sum_{i=1}^{N}\left(x_{i} - \bar{x}\right)^{2}
\end{align}

で与えられますが,この推定量は一致推定量(標本が十分大きくなると,真値である {\sigma^{2}}に収束する)であることも知られています.ここで, {\bar{x} = \sum_{i=1}^{N} x_{i}/N}であり,標本平均です.

推定量は確率変数の関数であることから,確率変数がばらつきを持つ以上,推定量もばらつきを持ちます.しかしながら,一致推定量は,標本サイズ {N}が大きくなると,真値からのばらつきはどんどん小さくなり,ばらつきが {0}に収束していくということを意味します.式で書くと,次のようになります.

任意の {\varepsilon > 0}に対して,次の式が成立する.


 P(|\hat\sigma^{2} - \sigma^{2}| > \varepsilon ) \rightarrow 0, \qquad n \rightarrow \infty

これは,「 {\hat\sigma^{2}} {\sigma^{2}}に確率収束する」ということを述べています.最後に,これを視覚的に表した図を示して終わろうと思います.以下のコードを実行します.これは,標本数がn.candの時のときの,標本分散をj.num(ここでは50回)計算して,それをプロットするということを行っています.これによって,標本数に対して推定量のばらつきがどのように変化するかが可視化されます.

set.seed(20160607)

true = 2
n.cand = c(2,3,5,10,20,50,100,200,500,1000,2000,10000)
j.num=50
res = matrix(0,length(n.cand),j.num)


for(i in 1:length(n.cand)){
	for(j in 1:j.num){
		n = n.cand[i]
		x = rnorm(n,mean = 0 ,sd = true)
		res[i,j] = var(x)
	}
}

for(j in 1:j.num){
	if(j==1){
		plot(log(n.cand),res[,j],type="l",ylim=c(min(res),max(res)),main="Convergence of Sample variance",ylab="estimated variance",xlab="log of sample size",xaxt="n")
		axis(1, at=log(n.cand), labels=n.cand)
	}else{
		points(log(n.cand),res[,j],type="l")
	}
}
abline(h=true**2,col=2)

結果の図

f:id:tomoshige_n:20160607231830p:plain
見ていただくとわかるように,徐々に赤線で書いた真値へ収束しているのが見て取れます.
やはり,視覚化するとわかりやすいです.内容はここまでですが,少々告知をさせてください.

宣伝

6月9日発売の岩波データサイエンスvol.3に「傾向スコアを用いたバント効果の解析」というテーマで寄稿させていただきました.大リーグのデータを用いてノーアウト1塁における犠牲バントの効果の検証を行った内容を掲載しております.本全体としましては(私も章立てしか存じませんが)調査観察データを用いた統計解析などに役立つ内容となっています.本は,基礎パートと,事例パートに分かれており,特に前半の星野先生の「統計的因果効果の基礎―傾向スコアと操作変数」と合わせて読んでいただければ,私の内容はわかりやすいかと思います.是非,手にとってご覧下さい.

公式:岩波データサイエンスVol.3
サポートページ:岩波データサイエンス

【2日目】統計を学ぶ人のための測度論(1週間限定独りリレーブログ)

こんにちは,2日目の記事はいろいろ悩みましたが,「統計のための測度論」ということで書いてみようかと.最初に断っておきますが,「理論的厳密さ」よりも,「直感的理解」を優先して書きますので,その辺り気持ち悪い人は,Wikipedia/数学書(最後の参考文献)などを参照ください.

さて,測度論といえば,Twitterをみている限り,勉強会で統計を勉強し始めた人が「本格的に避けたい」分野になっているような気がします.その実情が垣間見えるのは,こちら(※逆に,統計やってるのに測度知らないとか...みたいなことを書いてる人もいて,gkgkbrbrしました(´・ω・`)).

twitter.com

数学を専攻していた学部時代の僕でさえ正直なところ,統計やるんだからなんで必要なんだ?と思っていた時期があるぐらいですから,統計を知っておきたい/勉強を始めたい!という方に取って,これほど負担になっている分野はないことでしょう...

そこで,統計のための測度論というテーマでブログを書いてみます.できる限り,分かりやすく,説明する予定ですが,時間の関係上粗くなる箇所もあるかもしれませんのでご了承ください.

測度は「統計を勉強するとき」どこで使うのか?

勉強するには,モチベーションが大事なので,まずは「統計やってて,どこで測度ってでてくんねん?」を解消しましょう.誤解を恐れずにいえば,ビジネス書程度の統計の理解が必要な方(平均/分散/標準偏差/検定の基本事項などを頭に入れておくべき方)にとっては,この分野に時間を割くのは効率が良いとはいいがたいかなと思います.しかし,統計の世界をより詳しく,「確率って何だ?」「分布って?」「分布の期待値計算するとき,積分の中にf(x)とかp(x)とか書いてるけど,あれどっから来たし?」などなど疑問を晴らしながら進むためには,知っておいて損はありません.

皆さん大好き,MCMCだって,測度がわかっていなければ「アルゴリズムで得られるサンプルが,目的の事後分布からのサンプルになっている」なんて証明することすらできないんです(※Total Variation Normという,確率測度の差のノルムの上界が0に収束するということを示したりするのですが).

と,こんな感じで統計の確率分布などに触れて,細かいところが気になり始めると「測度」というのは切っても切り話せないお友達になるわけです.で,一口に測度論を勉強しましょうといいましても,数学的な基礎の基礎について復習をしていたら間に合いません.そこで,特に「お友達」になるべき「測度の基礎」と「確率測度」という測度に限ってここでは話していくことにします.

測度とは?

まず,測度とはその名の通り「何かを測る」ためのものだということです.
で,具体的に何を測るのかといえば,統計では特に「事象」の「起こりやすさ」を,[0,1]の値で表される「確率」に対応させて「測る」わけです.
(※単に測度といわれた場合には,[0,1]に対応させる必要はありません.)

そこで,統計て使われる「確率測度」の測度の感覚を掴むために,天気を例にとって考えます.

A を「雨が降る」という事象だと考えましょう.すると,この事象の起こりやすさをμという測度で測るというのは,「μ(A) = なんらかの値を対応づけること」と考えることができます.しかしながら,ここで問題になるのは「起こりやすさ」というのは「相対的な概念」であるという点です.つまり,「晴れる」「曇る」「雪が降る」などのそれぞれの事象があってはじめて,「雨が降る」という事象の「起こりやすさ」を定義することができます.

そこで,次のように考えれば良さそうです.

まず,「天気」のおこりうる全ての事象を含む集合を考えます(今回は簡単のため3つの事象で,晴れ,曇り,雨にしましょう).これを X = {晴れ,曇り,雨} とおいておくことにします.さらに仮定として,お天気は3つの事象のいずれか1つに定まるとして,かぶりがないとします.つまり,晴れでもあるし,曇りでもあるという状況がないと考えます(分かりにく場合には,1つのサイコロを振ったときが1の目と2の目は同時に出ないのと同じことです).

さて,このときまずはXそのものの集合を測ることを考えます.これは実際には,何にしても良いんです.{\mu(X) = 100}としてもいいし,{\mu(X) = 1}としても構いません.ただし,{\mu(X) = 0}や,{\mu(X) = -1}などとするのは,直感に反します(※理由は,後で分かりますが,起こりやすさの程度の議論をしていて,どれか一つが起こるはずのすべての事象の集合の大きさが0や負の数というのは,直感的におかしそうですよね).よって,我々としてはμ(X)に対して正の値を対応づけることにしましょう.ここでは,{μ(X) = a > 0}という値を対応づけます.

次に,{晴れ}という事象が,どの程度起こりやすいのかを考えます.ここで,μ({晴れ}) = b という値を対応づけます.ここでも,起こりやすさの議論をしているのでbは負の数にはなりません.一方で,b=0というのはあり得ます.つまり,{晴れ}になることが絶対になければb=0だからです.一方で,b > a となるのはおかしい気がします.つまり,天気をみたときに{晴れ,曇り,雨}のいずれかであるという「起こりやすさ」と,{晴れ}であるという起こりやすさは,明らかに現実的には前者の方が起こりやすいわけですから,b は aよりも小さいか,等しいでなければおかしいわけですね.

同様に考えると,{ \mu(\{曇り\}) = c,\mu(\{雨\}) = d},としておけば{b,c,d}はそれぞれ同じ性質を満たすはずです.では,最後に,先ほど「晴れ,曇り,雨」はどの2つも同時には起こらないといいました.つまり,{晴れ},{曇り},{雨}のどれか一つが起こる「起こりやすさ」の和,{晴れ, 曇り, 雨}の事象が起こる起こりやすさと等しいと考えるのが自然です.よって,

$$ a = μ(X) = μ(\{晴れ, 曇り, 雨\}) = μ(\{晴れ\}) + μ(\{曇り\}) + μ(\{雨\}) = b + c + d $$

となります.ここで,μ(X)にはどのような値を対応させても良かったわけですから,両辺をaで割った値を対応させることにすれば

$$1 = μ(X) = μ(\{晴れ, 曇り, 雨\}) = μ(\{晴れ\}) + μ(\{曇り\}) + μ(\{雨\}) = b/a + c/a + d/a$$

となります.このようにすれば「どれか一つが起こる」という起こりやすさは1であり,それぞれの起こりやすさは上記の足し算のそれぞれに対応します.ここで,{a > b,c,d}ですから,{ μ(\{晴れ\}),μ(\{曇り\}),μ(\{雨\})},はそれぞれ0以上1以下の値に対応づけられます.

これは,もしや,,,「確率」なんじゃないか!このように考えれば,確率と測度が結びつくわけですね.
では,もう少しこの議論を数学的に落とし込んでみましょう.

測度の定義を述べよう

上記の例は,あくまで「例」なわけで,実際に参考書においては「測度」を定義する際に,様々な言葉が使われています.ただ,我々の考えたい「測度」というのが,「ある事象」に対して「起こりやすさ」を,決めてやるものであるということは理解できると思います.これを数学的に落とし込んでみましょう.

ただ,上記では「起こりやすさ」と書いていますが,この「測度」の考え方はもともと,面積、体積、個数といった「大きさ」に関する概念を精緻化・一般化したものであることは指摘しておく必要があります.ここでいう「起こりやすさ」というのも,「ある事象」の「起こりやすさ」という「大きさ」の1つであると考えているわけですね.

さて,Wikipediaによると,測度とは以下のように定義されます.

集合 {X} の部分集合からなる完全加法族 {A} 上で定義される可算加法的測度 {\mu} とは拡張された区間 {[0, \infty]} に値を持つ(つまり、無限大も許す非負値の)関数であって、次の性質を満たすもののことである:

1.空集合の測度は 0 である。
$$ \mu(\emptyset) = 0. $$

2.完全加法性(可算加法性): {E_1, E_{2}, E_{3}, \cdots} がどの二つも互いに共通部分を持たない {A} に属する集合の列ならば
$$\mu\left(\bigcup_i E_{i} \right) = \sum_{i} \mu(E_{i}) $$
A の元は可測集合 (measurable sets ) と呼ばれる。 また、 数学的構造 {(X, A, \mu )} は 測度空間 (measure space ) と呼ばれる。 次の性質は、上の定義から導かれるものである:

{X}は「標本空間(sample space)」と呼ばれます.これは,ある実験(試行)を行ったときに, 起こり得る全ての結果の集合を指しています.上記の天気であれば,明日の天気を予想しましょう!という実験を行ったときに起こりうる全ての結果の集合であり,{X = \{晴れ, 曇り, 雨\}}のことになります.この結果の集合は「有限個」であっても,なくても構わないです.

次に,完全加法族{A}です(可算加法族/シグマ加法族と呼ばれることもあります).言葉としては難しいですが,起こりうる事象の冪集合{2^{X}}(の正確には,部分集合)で表現されるものであると考えてください.冪集合{2^{X}}といわれた場合には,これは{X}の任意の部分集合を元とする集合の族を指しています.これはつまり,{X = \{晴れ, 曇り, 雨\}}に対しては,次のように求めることができます.

$$ 2^{X} = \{ \{\emptyset\}, \{晴れ\}, \{曇り\}, \{雨\}, \{晴れ, 曇り\}, \{晴れ, 雨\}, \{曇り, 雨\}, \{晴れ, 曇り, 雨\}\} $$

ここで,どうして「完全加法族」というのが必要なの?という疑問が出てきます.標本空間だけあれば十分な気がするけど?と思うのですが...まずは,上記で定義される「完全加法族」というのがどのような性質を持つのかを確認します.

集合 {X} とその上の冪集合 {2^{X}} に対し、{X} の部分集合族 { \Sigma \subset 2^{X} }{X} 上の完全加法族であるとは、以下3つの性質を満たすものである

{\Sigma} は空でない: 少なくとも一つの {A \subset X}{\Sigma} に属する。

{\Sigma} は補演算に関して閉じている: {A}{\Sigma} に属するならば、その補集合 {X\A}{\Sigma} に属する。

{\Sigma} は可算合併に関して閉じている: {A_{1},A_{2}, A_{3}, \cdots}{\Sigma} に属する集合の列ならば、それらの合併 {A = A_{1} + A_{2} + \cdots}{\Sigma} に属する。

一方で,測度というのは「関数」なわけですから,どんな「範囲」で定義される関数なのかがとても大事になります.ここで,測度というものを完全加法族上の関数として考えてやることで直感に反しないように定義することができます.先ほどの議論からもわかるように,測度というのは,満たしてほしい重要な性質があります.

それは,まず空集合{\emptyset}に対しては{0}となることです.確率でいえば,ある試行を行った際,その結果が空集合に属するというのは「起こることはない」わけですから,空集合に対しては{0}であるように定義します.これが,上記の定義の1つ目です.

1.空集合の測度は 0 である。
$$ \mu(\emptyset) = 0 $$

次に,特にさきほど確率を考えているときには,{晴れ, 曇り}という事象が起こる「起こりやすさ」である{\mu(\{晴れ, 曇り\})}というのは,{晴れ}と{曇り}は同時に起こらない場合には,その「起こりやすさ」は$\mu$({晴れ}) + $\mu$({曇り})と分解しても,同じになってくれないと困るわけですね.そして,ここで先ほどの可算加法性が大事になるわけですが,{晴れ}と{曇り}という2つの事象の上で関数$\mu$が定義されていても,$\mu$が {晴れ, 曇り}上で定義されていないと困ります.そういった意味で,測度を可算加法族上での関数と考えるのは,目的にあっているような気がします.そして,確率を考える上で何より大事なのは,同時に起こらない2つの事象の和集合は,2つの事象の和に分解できるということです.$\mu$({晴れ, 曇り}) = $\mu$({晴れ}) + $\mu$({曇り}).これを,一般的に書いたものが以下の性質ということになります.

2.完全加法性(可算加法性): {E_1, E_{2}, E_{3}, \cdots} がどの二つも互いに共通部分を持たない {A} に属する集合の列ならば
$$\mu\left(\bigcup_i E_{i} \right) = \sum_{i} \mu(E_{i}) $$

これが「測度」というものです.最後に,用いたい「確率」という概念に対応させるためには,次の条件を追加してやれば良いわけですね.

全空間{X}の測度は,{\mu(X) = 1}である.

これで,確率測度が得られたことになるわけです.一般的に確率測度といわれた場合にはこのような{\mu}{P}として表すことが多く,このように記述されれば,「あ,確率のこといってるんだ!」とわかりますよね.

まとめー測度論とは結局「大きさ」を「測る」こと

結局,測度とは,「大きさ」を「測る」ことを意味しています.その大きさに「起こりやすさ」という意味を持たせると,そこに「確率」という概念が登場してくるわけです.統計の下地にある「測度」という概念は,ともすれば見落としがちですが,統計を学んでいくと「測度」について考えざるを得ないという状況はよく遭遇しますし,「測度」という概念を知っているだけで,読める本や論文の幅も広がります.上記の例では「離散」かつ「有限」な事象の可算加法族上の測度について話していますが,実際には「連続空間」上の「無限」な事象の可算加法族上でも測度は定義されます.そして,統計をやっている人が,きっと最初に学習する「正規分布」というのは,その測度の1つであり,何気なく「正規分布の仮定」でいいや!と言っているのは,その事象の起こり方について「自分が知っている」と仮定したもとで,事象の「起こりやすさの尺度」である「測度を定義している」ということに対応しています.また,目的変数が離散0か1なら「ロジスティック回帰」というのも,結果変数に対して「ベルヌーイ分布」を仮定していて,同じように起こりやすさの尺度を自分たちで無意識に仮定しているんですね.

と,こんな感じで,測度という概念はいろんなところに登場しておりまして,普段統計を使うときには気にしないでもまぁ大丈夫なものだと僕は思っておりますが,緻密に議論をしたい,統計をやっていて「ここは何だ?」と思い始めてしまった,そんなときには是非「測度」に立ち返ってみようかな?と思ってもらえれば嬉しいかなと思いました.

僕も,最初はきらいだった測度ですが,最近は友達ぐらいにはなれたかな?と思っています.周囲に,こういう数学の概念で分からないところがあれば尋ねれば「理解して,実際に使っている人がいる」というのが,大学院で統計や数学を専攻している1つのメリットかも,,,と改めて思った次第です.

参考

この本とってもいい本です.最近でも見直すレベルです.ちょっと高いですが苦笑

  • Real Analysis

Amazon.co.jp: Real Analysis: Modern Techniques and Their Applications (Pure and Applied Mathematics: A Wiley Series of Texts, Monographs and Tracts) 電子書籍: Gerald B. Folland: 洋書


これは,友達が薦めてくれた本で,僕が測度論をわかった!と思わせてくれた1冊.

  • Real and Complex Analysis

Amazon.co.jp: Real & Complex Analysis: Rudin: 洋書


Wikipediaは偉大です.

測度論 - Wikipedia
冪集合 - Wikipedia
完全加法族 - Wikipedia

p.s.

今年も,アドベントカレンダーが盛り上がっています.特に,leafletは使ってみたいなと思っていたので参考にして何かやってみるつもりです.
qiita.com

次回は

確率密度関数/確率関数,分布関数について書いてみます.明日は,おやすみします.

【1日目】1週間限定独りリレーブログ

今年も,もうすぐ終わりですね.そして,更新しなきゃ!と思って放置していたこのブログも遂に放置して7ヶ月が経とうとしています.真面目に更新しろよ,自分.ということで,年末こそは更新するぞを息巻いて,1日目の更新をします(※このブログも,現在月に1万Viewぐらいは頂いているようで,嬉しい限りです).毎年,この時期になると,「指が冷たいから,数式書くのやだ...」とか,「やばい,論文が間に合わない」と焦っているのですが,今年もそんな感じでバタバタです.

2015年の研究活動の振り返り

2015年は,どんなことをしていたのかというと主に対外発表をさせて頂きました.
・初国際学会で発表(8月 シアトル)
・初国内学会で発表(9月 岡山)
・初研究集会で発表(11月 東京)

特に思い出深かった8月の国際学会の発表を書いてみようかなと.

国際学会(8月)

JSM2015(Joint Statistical Meeting 2015)に行かせて頂いて,発表をさせて頂きました.
www.amstat.org

これが,自分にとっては2015年一番印象に残っている出来事かなと思います.参加者が7000+人,セッション数700というみたこともない規模に圧倒されながら,そして自分のセッションでは足が震えながら発表しました(英語の添削などもいろんな方にお手伝いいただきました.日本語のニュアンスを,英語に翻訳するのは難しいことを実感しましたし,もっと〆切前に頑張らねばというコトを痛感した次第.).

(初参加の人向けのセッション写真)
f:id:tomoshige_n:20150809124554j:plain

(エントランス周辺の人だかり)
f:id:tomoshige_n:20150810122939j:plain

僕の発表は,タイトルの通り微小粒子状物質に関連するお話で,現状の解析手法で懸念されることや,空間上のPM2.5の分布のBayes予測の話をしました.

また,JSM2015では僕が意識していたからかもしれませんが,タイトルだけみても,データの「Dependence」に関連する発表は比較的多かったように思います.解析では,よく置かれる仮定として,それぞれのサンプルは「独立同一に...」というようなものがありますが,実際に「遺伝子」や「リスク管理」などの分野では,サンプル間に相関/関連性が認められることが多く,このような場合には,どのような解析手法を用いるのが適切なのか?というのがいくつかのセッションで発表/議論されていました.

発表の分野としては,「機械学習」もありましたし,「経営」「マーケティング」「因果推論」「スポーツ(メジャーリーグの投手の解析)」などもあったようです.同時進行でセッションは進んでいるため,全てに出ることは叶わないのですが,アブストラクトをみていても,これはおもしろいな!と思えるものがいくつもありました(終了後にも,みれるようにしてほしいのですけどね笑).

あと,EXPOというのも開催されていて,企業ブースがあったりします.SAS, R-Studio, Revolution R, Springer, CRCなどなど,統計やってる人に取ってはお世話になってる企業がずらりと揃っていました.ステッカーいろいろ頂きまして,Macbook Airに貼らせてもらっています♪

f:id:tomoshige_n:20150810144101j:plain

p.s.

セッション終了後には,近くのマーケットに御飯を食べにいったりしましたが,スタバの1号店が近くにありました.
あと,アメリカで食べるSubwayはボリュームが大きくて,格別でした!(/∇\*)
f:id:tomoshige_n:20150811184630j:plain

研究以外のこと

Stanを使い始めた.

そういえば,Stanを本格的に使い始めました.なんというか,自分でMCMC書くのがちょっと面倒になってきまして,モデルをいくつも試すなら,やっぱりStanかなぁと思って,Stanを岩波データサイエンスとStan Manual v2.8.0を片手にいろいろ書いてみたのですが,確かにうまくサンプルとれてそうだけど書き方が下手なのか,とにかくサンプリングが遅い気がして,悶々としています.ただ,とっても楽です♪.あと,使えることは分かったけど,これ理論的な証明どうなってるんだろ?というところで悩み始めておりまして,論文を読みあさっております.
岩波データサイエンス

今年の猛省.

これは,「オンラインゲーム」(オルクスオンライン)を様々な言い訳をつくってやり過ぎたという一言につきますね.とりあえず,ゲームはリタイアして,研究します(・ε・).
aurcus.jp

次回

明日は,年末年始にできる簡単な統計の復習を書いてみます〜!

パッケージを使わないで、一般化線形混合モデルのMCMCアルゴリズムを1から作る.

こんにちは.
金曜日の夜になり、激しめの睡魔に襲われております.

先日のこちらの記事で公開したスライドの後半にあるシミュレーションで、地域差を考慮したPoisson - Normalモデルを構築しているのですが、そのコードを載せておきます。tomoshige-n.hatenablog.com

MCMCを自分で設計

最近だと、MCMCはStanだったり、BUGSだったり、色々と便利なツールもできてきていて、パッケージを用いて推定している人も多いみたいです。が、僕は、なんというか、中で何してるのかが気になって、提案分布の設計とか実際のところ、どうしてるのかとか、ステップ幅をどうやって決めてるのかとか、自動的にやられてしまうのはどうも気持ち悪いので、自分で書いてみましょうということになってしまいがちなんですよね。

ということで、ポアソン分布で、個体差を考慮できるようなMCMCアルゴリズムを書いてみました.ちなみに、大変申し訳ないのですが、こちらの記事のMCMCのコードには、いくつか間違いがありますので(そのまま載せてはおきますが)、、ご注意ください.

tomoshige-n.hatenablog.com

書いたRのコード

まぁ、提案分布のステップ幅などは、きちんと調節しなくては行けないですし、収束の判定をするようには作っていないので、、、申し訳ないですが、その辺が必要な方はコードを書き加えていただければ...(追記)「パッケージ使わないで、、、」とかいいつつ、truncated normalとかはパッケージ使って計算してます.許してください.

#poisson regression

set.seed(20150513) #wed seminar date
library(dummies)
library(truncnorm)
library(MASS)
#simulation
a = 400
b = 400
c = 400
x = data.frame(
	x1 = c(runif(a,-2,2), runif(b,-1.8,2.2),runif(c,-3,2.4)),
	#x2 = c(rnorm(a,-0.5,1), rnorm(b,0,1),rnorm(c,0,1))
	x2 = c(rep(1,a), rep(0,b), rep(0,c))
)
x = as.matrix(x)
z_fac = c(rep("A",a),rep("B",b),rep("C",c))
z_fac1 = dummy(as.factor(z_fac))
z = z_fac1[,-1]
beta = c(0.8,-0.5)
u = c(rnorm(a,0,0.2),rnorm(b,0,0.3),rnorm(c,0,0.3))
lambda = exp(x%*%beta + u) #z%*%u は正規分布からの乱数と考えればよい。
y = rpois(a+b+c,lambda)

par(mfcol=c(1,2))
plot(x[,1],y,xlab="",ylab="",col=as.factor(z_fac))
plot(x[,2],y,xlab="",ylab="",col=as.factor(z_fac))

#個体差なし
y_norandom = rpois(a+b+c,exp(x%*%beta))
par(mfcol=c(1,2))
plot(x[,1],y_norandom,xlab="",ylab="",col=as.factor(z_fac))
plot(x[,2],y_norandom,xlab="",ylab="",col=as.factor(z_fac))


############################
#mcmc calculation from here#
############################

#1.proposal distribution

## regression coefficient parameter
betaprop = function(beta_before,step_width,step_width2,sigma,p,q){
	n = length(beta_before)
	#mu = c(beta_before[seq(1,p,1)],rep(0,q))
	mu = beta_before
	Cov = diag(c(rep(step_width,p),rep(step_width2,q)))
	sample = mvrnorm(1, mu = mu, Sigma = Cov)
	return(sample)
}

## area dispersion parameter
sigmaprop = function(sigma,sd_trunc,q){
	mu = sigma
	sample = rep(0,q)
	for(i in 1:q){
		sample[i] = rtruncnorm(1, a=0.05, b=Inf, mean = mu[i], sd = sd_trunc)
	}
	return(sample)
}

## to calculate proposal distribution probability (q(x_{n} | x_{n-1}))
prop_lograte = function(beta_after,beta_before,sigma_after,sigma_before,p,q,sd_trunc,step_width2){
	#beta_log = sum( - (step_width2/2)*((beta_after[seq(p+1,p+q,1)])**2))
	sigma_log = sum(
		-(1/2)*((sigma_after-sigma_before)/sd_trunc)**2 - log(1 - pnorm((0.05-sigma_before)/sd_trunc))
	)
	return(sigma_log)
}

#2. posterior

prior_logsigma = function(sigma,theta,k){
	result = sum((k-1)*log(sigma) - (1/theta)*(sigma))
	return(result)
}

#mixed effectの所のlogは計算する必要がある。
prior_logmixbeta = function(beta,sigma){
	#betaとsigmaは次元分あるので、最後にsumを取る
	result = sum(-(1/2)*((beta**2)*sigma) + (1/2)*log(sigma))
	return(result)
}

#improperな事前分布は収束しないこともあるので、betaのpriorを分散1000の正規分布にしてみた.
prior_beta = function(beta){
	#betaとsigmaは次元分あるので、最後にsumを取る
	result = sum(-(1/(2*1000))*((beta**2)) - (1/2)*log(10))
	return(result)
}

#正規化定数なしの事後分布の確率のlog
log_post = function(y,x,z,beta,sigma,theta,k){
	p = dim(x)[2]
	q = dim(z)[2]
	beta1 = beta[seq(1,p,1)]
	beta2 = beta[seq(p+1,p+q,1)]
	log_lambda = x%*%beta1 + z%*%beta2
	likeli = sum(log_lambda*y) - sum(exp(log_lambda))
	#prior calculation
	#beta given sigma
	beta_given_sigma_prior = prior_logmixbeta(beta2,sigma)
	#sigma
	prior_logsigma = prior_logsigma(sigma,theta,k)
	#beta
	#improperな事前分布を避ける
	beta_prior = prior_beta(beta1)
	result = likeli + beta_given_sigma_prior + prior_logsigma + beta_prior
	return(result)
}

#アルゴリズム
mcmc_pois = function(y,x,z,intercept = F,iter=1000,burn_in = 500,step_width=0.01,step_width2=0.01,theta=0.01,sd_trunc=0.01,k=0.01){
	p = dim(x)[2]
	q = dim(z)[2]
	initial_beta = c(rnorm(p),rep(0,q))
	initial_sigma = runif(q,0.5,1)
	beta.r = matrix(0,nrow=p+q,ncol=iter)
	sigma.r = matrix(0,nrow=q,ncol=iter) 
	log_likeli = rep(0,iter)
	accept = 0
	#初期値の設定
	beta.r[,1] = initial_beta
	sigma.r[,1] = initial_sigma
	log_likeli[1] = log_post(y,x,z,initial_beta,initial_sigma,theta,k)
	accept.r = rep(0,iter)
	alpha.r = rep(0,iter)
	for(i in 2:iter){
		#sigmaの候補を生成
		candidate_sigma = sigmaprop(sigma.r[,i-1],sd_trunc,q)
		#betaの候補を生成
		candidate_beta = betaprop(beta.r[,i-1],step_width,step_width2,candidate_sigma,p,q)
		#この値に基づくlog-postを計算する
		log_likeli_candidate = log_post(y,x,z,candidate_beta,candidate_sigma,theta,k)
		#棄却率alphaの計算
		alpha = exp(
			log_likeli_candidate - 
			log_likeli[i-1] - 
			prop_lograte(candidate_beta,beta.r[,i-1],candidate_sigma,sigma.r[,i-1],p,q,sd_trunc,step_width2) +
			prop_lograte(beta.r[,i-1],candidate_beta,sigma.r[,i-1],candidate_sigma,p,q,sd_trunc,step_width2)
		)
		#alphaに基づく条件分岐
		alpha.r[i] = alpha
		if(runif(1) < alpha){
			beta.r[,i] = candidate_beta
			sigma.r[,i] = candidate_sigma
			log_likeli[i] = log_likeli_candidate
			accept.r[i] = 1
		}
		else{
			beta.r[,i] = beta.r[,i-1]
			sigma.r[,i] = sigma.r[,i-1]
			log_likeli[i] = log_likeli[i-1]
			accept.r[i] = 0
		}
	}
	print_obj = sum(accept.r)*100/iter 
	print(paste("accept rate is ",print_obj,"%",sep=""))
	return(list(
		sample_beta = beta.r,
		burn_in_beta = beta.r[,-c(1:burn_in)],
		sample_sigma = sigma.r,
		burn_in_sigma = sigma.r[,-c(1:burn_in)],
		accept=accept.r,
		alpha = alpha.r,
		ll = log_likeli))
}

result = mcmc_pois(y,x,z,intercept = F,iter=1000000,burn_in = 200000,step_width=0.001,step_width2=0.001,k=0.01,theta=(1/0.01),sd_trunc=0.7)

あとは、自己相関とか、収束先の分布などをプロットするために、以下のコードを使えばいいと思います。

#サンプルの自己相関の確認
par(mfcol=c(2,1))
acf(result$burn_in_beta[1,], type = "correlation", plot = TRUE, lag.max=400)
acf(result$burn_in_beta[2,], type = "correlation", plot = TRUE, lag.max=400)
par(mfcol=c(2,1))
acf(result$burn_in_beta[3,], type = "correlation", plot = TRUE, lag.max=400)
acf(result$burn_in_beta[4,], type = "correlation", plot = TRUE, lag.max=400)
par(mfcol=c(2,1))
acf(1/result$burn_in_sigma[1,], type = "correlation", plot = TRUE, lag.max=800)
acf(1/result$burn_in_sigma[2,], type = "correlation", plot = TRUE, lag.max=800)

#自己相関が小さくなる適切な間隔毎にデータを引っこ抜く
sample_thin = function(x,lag){
	n = dim(x)[2]
	num = seq(1,n,lag)
	return(x[,num])
}

#上記の操作で得られるMCMCのサンプル
mcmc_sample_beta = sample_thin(result$burn_in_beta,400)
mcmc_sample_sigma = sample_thin(result$burn_in_sigma,400)

#事後分布
par(mfrow=c(2,2))
hist(mcmc_sample_beta[1,],xlab="",ylab="",main="beta1",breaks=15) #0.8
hist(mcmc_sample_beta[2,],xlab="",ylab="",main="beta2",breaks=15) #-0.5
hist(sqrt(1/mcmc_sample_sigma[1,]),breaks=50,main="sigmaB") #発生元は0.3
hist(sqrt(1/mcmc_sample_sigma[2,]),breaks=50,main="sigmaC") #発生元は0.3

あと、やはり予測値というか、値の信頼区間内にデータがどの程度入るのかが気になるので、予測値の90%信頼区間書いてみます。ここは、もっと効率よくできるはずですが、、、やる気スイッチが切れたのでこれで勘弁してください。

#上のMCMCサンプルデータを一つのマトリックスにまとめる
prdt_para = data.frame(
	beta1 = mcmc_sample_beta[1,],
	beta2 = mcmc_sample_beta[2,],
	sigma1 = 1/mcmc_sample_sigma[1,],
	sigma2 = 1/mcmc_sample_sigma[2,]
)

#推定値を保存するところ
prdt_para = as.matrix(prdt_para)
p = dim(x)[2]
q = dim(z)[2]
n = length(y)
r = dim(prdt_para)[1]
result2 = matrix(0,n,r)
result3 = matrix(0,n,r)

#各点の結果を予測する
for(i in 1:n){
	result_ind = rep(0,r) #各地点の結果を保存する
	for(j in 1:r){
		fix_eff = prdt_para[j,seq(1,p,1)]
		#rand_eff = mvrnorm(1,rep(0,q),diag(prdt_para[j,seq(p+1,p+q,1)],q))
		rand_eff = mvrnorm(1,rep(0,q),diag(c(0.26,0.28),q))
		mu = exp(t(x[i,]) %*% fix_eff + t(z[i,]) %*% rand_eff)
		result2[i,j] = mu
		result3[i,j] = rpois(1,mu)
	}	
}

#上側97.5%点
up95 = function(x){
	quantile(x,0.975)
}

#下側97.5%点
bo05 = function(x){
	quantile(x,0.025)
}

#ここから、図を描く作業
upper_line = apply(result3,1,up95)
bottom_line = apply(result3,1,bo05)
predicted_mu = apply(result2,1,mean)
mu_ord = order(predicted_mu)

plot(predicted_mu,y,col = as.factor(z_fac),pch=16,ylim=c(0,25))
abline(a=0,b=1)
points(predicted_mu[mu_ord],bottom_line[mu_ord],type="l",lwd=0.7)
points(predicted_mu[mu_ord],upper_line[mu_ord],type="l",lwd=0.7)
#upper line
loess_up = loess(upper_line[mu_ord] ~ predicted_mu[mu_ord])
xl = seq(0,20,0.01)
lines(xl, predict(loess_up,xl), col='red', lwd=2)
#bottom line
loess_bo = loess(bottom_line[mu_ord] ~ predicted_mu[mu_ord])
xl = seq(0,20,0.01)
lines(xl, predict(loess_bo,xl), col='red', lwd=2)

ちなみに、信頼区間的なものを書いておりますが、ポアソン分布は「左右非対称」なので、Highest posterior region(もっとも事後分布の値が大きい95%の区間)を求めるのは、難しいので、上記のような処理に落ち着けました...

結果でこんな図が書けます.

事後分布

f:id:tomoshige_n:20150515194554p:plain

信頼区間とか

f:id:tomoshige_n:20150515195036p:plain

まとめ

意外に自分でMCMCを書くのは楽しかったです.特に、計算を速くさせるための工夫とか、やっぱり大事だなと。特に、積はlogとれば和の計算というのは、数値計算では有効だなーと思いました。特にガンマ関数、階乗の計算はさせるべきじゃありません...MCMCのRのコードを書きながら、昔、MCMCについて素人の頃に、全然MCMCのコードを書いてる記事が見つからなくて困ったのを想い出しました。そんな人の少しでも助けになればと思いながら、まとめとさせてもらいます。

おすすめ文献

参考にした文献が今回はないので、読んでいる本を載せておきます.僕の研究分野の本です.

Hierarchical Modeling and Analysis for Spatial Data, Second Edition (Chapman & Hall/CRC Monographs on Statistics & Applied Probability)

Hierarchical Modeling and Analysis for Spatial Data, Second Edition (Chapman & Hall/CRC Monographs on Statistics & Applied Probability)

ベイズ推論とシミュレーション法の基本について

ベイズ推論とMCMC

今日セミナーで話をするスライドです.
基本的なベイズ推論と、シミュレーション法の枠組みを
自分なりにまとめたスライドになっています。

ベイズ推論は便利ですが、なかなか背景の理解が難しいので、
ベイジアンとノンベイジアンな手法の比較をしながら、
つたないながら説明を試みてみました.

参考になれば、幸いです.

www.slideshare.net