stanで平均値の差の推定。収束診断とMAP推定値の計算。LOOとWAICの算出。事後分布の可視化まで。

author: Unadon (見習い飯炊き兵) 動作環境:Mac OS Sierra 10.12.1; R version3.3.1; rstan 2.10.1

 

今回は、これができるようになります: 二群の平均値の差の推定

はじめに

rstanで色々ベイズしてみたい!という方に向けて、基本的な解析とその周辺Tipsをご紹介していきます。

今回は「平均値の差の推定」です。

とりあえずやってみて、やりながらわかっていく。

わかるからできるのではなく、「できるからわかる」方向での学びに貢献できたらと思います。

こまったら以下を御覧ください。だいたいここに行き着きます。

http://www.slideshare.net/simizu706/stan-62042940

目的

  1. Rとrstanで平均値の差の推定ができるようになる

  2. 収束診断、MAP推定値など、推定に関する必要な情報を出力できるようになる

  3. 結果の図を作成する。

ライブラリの読み込み

##################################################
#Stan : Gussian Difference
##################################################
invisible({rm(list=ls());gc();gc()})

#library and directory
library(rstan) 
library(ggplot2)
library(reshape2) #データ成形ツール
library(shinystan) #収束診断を含む色々な情報をまとめて知ることができる
library(loo) #Leave-One-Out クロスバリデーション、WAIC情報量基準の算出
library(tidyverse)  #データ成形最強ツール
library(plotly) #触れるグラフ
library(Rmisc)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

#
# 目的: 正規分布でデータの平均と標準偏差を推定し、
#         2群間の平均値の差を推定する
#


#Options-----------------------------------------------------
# 1. WAIC&LOO: 情報量規準WAIC&LOOの算出
# 2. Convergence Diagnosis: shinystanで収束診断 
# 3. Parameters Visualization of mu, sigma, Difference
# 4. MAP estimate
#------------------------------------------------------------

推定の設定とデータの読み込み

みなさまにいじってもらうのはここでしょうか。サンプリングの回数などを調整します。

データは今回irisを使いますが、ご自身のデータをcsvファイルから読み込んで頂いても構いません。 サンプリングの条件とデータのcsv読み込み、加えて変数名(データの列名)さえ指定すればOKです。

#Initaial Settings-------------------------------------------
#サンプリング数
iterations =c(2000)
#Warm-up数
warmup=c(1000)
#チェーン数
chains=c(3)
#thinning
thin =c(1)

#推定に使用するパラメータ,muとsigmaはベクトル
parameters=c("mu","sigma","Difference","log_lik") #log_likは必ず最後

#あとで結果を確認したいパラメータ
parameters2<-c("mu","sigma","Difference")

#データセットの読み込み(行は参加者,列は変数。Headerあり)
#dat<-read.csv("iris.csv")
#read data サンプルデータ
dat<-iris


#平均と標準偏差を求めたいデータの列名
varNames<-c("Sepal.Length","Petal.Length")

#列名でデータを抜き出してマトリックスにする
Dat<-data.matrix(dat[,varNames])
#------------------------------------------------------------

#Number of variables and Subjects
Nvar<-length(Dat[1,])
Nsub<-length(Dat[,1])

#data and parameters to be passed on stan
datastan=list(Nvar=Nvar,Nsub=Nsub,Dat=Dat)


使うデータは2列の150行。表形式です。

>head(Dat)

     Sepal.Length Petal.Length
[1,]          5.1          1.4
[2,]          4.9          1.4
[3,]          4.7          1.3
[4,]          4.6          1.5
[5,]          5.0          1.4
[6,]          5.4          1.7

Stanのモデルコード

以下はStanのモデルコードです。

Sepal.LengthとPetal.Lengthのデータをマトリックスとして扱い、それぞれの平均値を求めています。 途中、Differenceを計算しているのは、毎回のサンプルにおけるSepalとPetalの差です。この差を集めて、差の事後分布を得ていきます。


GaussianDifference<-"
data {
    int <lower=0>Nvar;
    int <lower=0>Nsub;
    matrix [Nsub,Nvar] Dat;
}

parameters {
    vector [Nvar]mu;
    vector <lower=0>[Nvar]sigma;
}

transformed parameters {
    real Difference;
 Difference = mu[1]-mu[2];
}

model {
    for (i in 1: Nsub){
        for (j in 1:Nvar){
            Dat[i,j] ~ normal(mu[j],sigma[j]);
        }
    }
}

generated quantities {
    matrix [Nsub,Nvar]log_lik;
    for (i in 1:Nsub){
        for (j in 1:Nvar){
            log_lik[i,j] = normal_lpdf(Dat[i,j]|mu[j],sigma[j]);
        }
    }
}

"


サンプリングの実行と推定結果の要約

サンプリングを実行します。結果の出力は、Summaryのところのコマンドです

#StanFit
fit<-stan(model_code=GaussianDifference,data=datastan,pars=parameters,
              iter = iterations,chains = chains,thin = thin)

#Summary
Summary<-summary(fit)$summary[c(1:c(Nvar*Nvar+1)),c(1,3,4,8,9,10)]

推定されました。Differenceの95%確信(信用)区間が0を超えているので、差があると言えそうですね。

print(Summary)

              mean         sd      2.5%     97.5%    n_eff      Rhat
mu[1]      5.8422909 0.07048578 5.6997652 5.9808916 2597.486 0.9999456
mu[2]      3.7633578 0.14844683 3.4707780 4.0540803 2739.429 0.9997349
sigma[1]   0.8363946 0.05077964 0.7456214 0.9474945 2876.424 1.0000821
sigma[2]   1.7831513 0.10920370 1.5903404 2.0028148 2207.373 1.0003378
Difference 2.0789331 0.16138221 1.7671881 2.4055275 2647.713 0.9997604

収束診断

まずやるべきは「この推定が上手く行ったのか」確認する作業です。 Gelman-Rubin testを自動で実施してくれるshinystanを起動してみましょう。

launch_shinystan(fit)

次のような画面がでるので、DIAGNOSEをポチ。

shinystan.png

RhatやMCSEといった指標のところをみて、”None”が表示されていたらテストはクリア、推定は上手くいっています。

shinystanDiagnosis.png

結果の取り出しと可視化

shinystanをいじってもらうと、推定結果の事後分布やサンプリング過程など、色々必要な情報が得られます。

それでも良いのですが、軸を変えたり、色を変えたりなど、図の見せ方を変えたいときは、自分で作図する必要があります。

やってみましょう。 ややこしいコードに見えるかもしれませんが、要点を押さえるとだんだんわかってきます。

何があればベイズの結果の図が書けるのか?

  1. パラメータの名前
  2. 対応したMCMCサンプル

この2つだけです。

あとは,ggplotに求められる整然データ(tidy data)になるようなデータの成形と、 ggplotの書き方の問題です。



#パラメータのggplot----------------------------------------

#Extract(全サンプルの取り出し)
d.ext<-rstan::extract(fit)    
#上で決めていた必要なパラメータの名前をつかって、当該パラメータのサンプルだけを抜き出す
df.ext<-data.frame(d.ext[c(parameters2)]) 

#melting to tidy data ここは気にしない。
df.long<-melt(df.ext,idvar=c())

#col name separation ここも気にしない
df.comp<-df.long %>% separate(variable, c('parameter', 'varNames'))

#replace Num of parameter to varNames 無視無視
for(i in 1: Nvar){
  df.comp$varNames<-replace(df.comp$varNames,df.comp$varNames==as.character(i), varNames[i]) 
}
df.comp$varNames<-replace(df.comp$varNames, which(is.na(df.comp$varNames)), "Difference")  
df.comp$varNames<-as.factor(df.comp$varNames)

#ここまでで作図に必要なデータの加工が終了


####ここから図を作っていく

#theme set 真っ白なキャンバスを用意
ggplot()+theme_set(theme_bw(base_size=18,base_family="HiraKakuProN-W3"))

#色々書き込んでいく

#データはdf.compをつかって、x軸は変数名(列名)value, y軸は確率密度で、線の色は変数名で変えて、塗りつぶす色も変数名で変えることを指定
g <- ggplot(df.comp,aes(x = value,y = ..density..,colour=varNames,fill=varNames))

#ヒストグラムをキャンバスにかきこむ。αは透過性
g <- g + geom_histogram(position = "identity",alpha = 0.5)

#密度曲線をキャンバスに重ね書き
g <- g + geom_density(stat = "density",position = "identity",alpha = 0.4)

#facetする。キャンバスをパラメータの種類で分ける
g<-g+facet_wrap(~parameter,scales="free")

#色わけは”Set2”で。
g<-g + scale_color_brewer(palette="Set2")
g<-g + scale_fill_brewer(palette="Set2")

#タイトルは
g<-g+ggtitle("Posterior Distributions")

#axis 軸のタイトルとテキストサイズをいじる
g <- g + theme(axis.title.x = element_text(size=14), #axis title letter size
               axis.title.y = element_text(size=14)) 
g <- g + theme(axis.text.x = element_text(size=8),  #axis text letter size 
               axis.text.y = element_text(size=8))   
g <- g + theme(legend.title = element_text(size=14), #legend size
               legend.text = element_text(size=10),  #legend text size 
               legend.key= element_rect(colour = "white")) # legend mark-box colour
#出力
g


#インタラクティブに表示したい場合はこちら
#if you need interactive plot, 
ggplotly(g)


で、よいしょっと。 最頻値あたりにカーソルをあわせて値を見ておいてください。

事後分布の最頻値(Maximum a Posteriori: MAP推定値)の計算

事後分布に対称性がなく、右に歪んでたりするばあいは、推定平均値ではなく、確率密度の最も高いMAP推定値を代表値として用いることがあります。 Rで計算します。


#MAP estimate-------------------------------------------------------
MAPs  = c(rep(NA,c((length(parameters2)-1)*(length(parameters2)-1)+1)))
for (i in 1:c((length(parameters2)-1)*(length(parameters2)-1)+1)){
    MAPs[i]    <- density(df.ext[,i])$x[which(density(df.ext[,i])$y == max(density(df.ext[,i])$y))]
}
names(MAPs) <- c(colnames(df.ext))

でてきました。どうやらPetalよりSepalの方が2くらい長いらしい。

print(MAPs)
      mu.1       mu.2    sigma.1    sigma.2 Difference 
 5.8309381  3.7810705  0.8308704  1.7394855  2.0690658 

情報量基準の算出

今回の平均値の差の推定では、まぁ不要だと思いますが、 ベイズでは、幾つかのモデルをデータに適用して、どれがいいのか(予測の良さや当てはまりの良さ)を比べることがあります。相対比較です。

今回は「予測のよさ」の指標であるLOOとWAICを算出します。 モデルコードの中で推定していた対数尤度log_likを使って計算します。

Petal.LengthとSepal.Lengthのそれぞれについて、 正規分布を適用したことによる予測的な当てはまりの良さを求めます。 要するに、指標が2個出てくるということです。 正規分布ではなく別の分布を当てはめたときの場合などと比較してください。

今回はそれを行いませんが。

#情報量規準WAICとLOOの算出
#WAIC and Leave one out cross validation--------------------------------------------
WAIC<-c()
for (i in 1:Nvar){
    res<-waic(extract_log_lik(fit)[,c(Nsub*(i-1)+1):c(i*Nsub)])
    WAIC[i]<-res$waic
}

#Leave-One-Out cross validation
LOO<-c()
for (i in 1:Nvar){
    res2<-loo(extract_log_lik(fit)[,c(Nsub*(i-1)+1):c(i*Nsub)])
    LOO[i]<-res2$looic
}

結果です。

> print(WAIC)
[1] 372.0823 598.7035

> print(LOO)
[1] 372.0820 598.7041

最後に

今回ご紹介した平均値の差の推定コードをまとめて掲載しておきます。

次回は線形モデルと事後予測チェックのコードをアップします。

では、Enjoy Stan!!

##################################################
#Stan : Gussian Difference
##################################################
invisible({rm(list=ls());gc();gc()})

#library and directory
library(rstan) 
library(ggplot2)
library(reshape2) #データ成形ツール
library(shinystan) #収束診断を含む色々な情報をまとめて知ることができる
library(loo) #Leave-One-Out クロスバリデーション、WAIC情報量基準の算出
library(tidyverse)  #データ成形最強ツール
library(plotly) #触れるグラフ
library(Rmisc)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

#
# 目的: 正規分布でデータの平均と標準偏差を推定し、
#         2群間の平均値の差を推定する
#


#Options-----------------------------------------------------
# 1. WAIC&LOO: 情報量規準WAIC&LOOの算出
# 2. Convergence Diagnosis: shinystanで収束診断 
# 3. Parameters Visualization of mu, sigma, Difference
# 4. MAP estimate
#------------------------------------------------------------



#Initaial Settings-------------------------------------------
#サンプリング数
iterations =c(2000)
#Warm-up数
warmup=c(1000)
#チェーン数
chains=c(3)
#thinning
thin =c(1)

#推定に使用するパラメータ,muとsigmaはベクトル
parameters=c("mu","sigma","Difference","log_lik") #log_likは必ず最後

#あとで結果を確認したいパラメータ
parameters2<-c("mu","sigma","Difference")

#データセットの読み込み(行は参加者,列は変数。Headerあり)
#dat<-read.csv("iris.csv")
#read data サンプルデータ
dat<-iris


#平均と標準偏差を求めたいデータの列名
varNames<-c("Sepal.Length","Petal.Length")

#Data shaping
Dat<-data.matrix(dat[,varNames])
#------------------------------------------------------------

#Number of variables and Subjects
Nvar<-length(Dat[1,])
Nsub<-length(Dat[,1])

#data and parameters to be passed on stan
datastan=list(Nvar=Nvar,Nsub=Nsub,Dat=Dat)


GaussianDifference<-"
data {
    int <lower=0>Nvar;
    int <lower=0>Nsub;
    matrix [Nsub,Nvar] Dat;
}

parameters {
    vector [Nvar]mu;
    vector <lower=0>[Nvar]sigma;
}

transformed parameters {
    real Difference;
 Difference = mu[1]-mu[2];
}

model {
    for (i in 1: Nsub){
        for (j in 1:Nvar){
            Dat[i,j] ~ normal(mu[j],sigma[j]);
        }
    }
}

generated quantities {
    matrix [Nsub,Nvar]log_lik;
    for (i in 1:Nsub){
        for (j in 1:Nvar){
            log_lik[i,j] = normal_lpdf(Dat[i,j]|mu[j],sigma[j]);
        }
    }
}

"


#StanFit
fit<-stan(model_code=GaussianDifference,data=datastan,pars=parameters,
              iter = iterations,chains = chains,thin = thin)

#Summary
Summary<-summary(fit)$summary[c(1:c(Nvar*Nvar+1)),c(1,3,4,8,9,10)]



#パラメータのggplot----------------------------------------

#Extract(サンプルの取り出し)
d.ext<-rstan::extract(fit)    
#上で決めていた必要なパラメータの名前をつかって、当該パラメータのサンプルだけを抜き出す
df.ext<-data.frame(d.ext[c(parameters2)]) 

#melting to tidy data ここは気にしない。
df.long<-melt(df.ext,idvar=c())

#col name separation ここも気にしない
df.comp<-df.long %>% separate(variable, c('parameter', 'varNames'))

#replace Num of parameter to varNames 無視無視
for(i in 1: Nvar){
  df.comp$varNames<-replace(df.comp$varNames,df.comp$varNames==as.character(i), varNames[i]) 
}
df.comp$varNames<-replace(df.comp$varNames, which(is.na(df.comp$varNames)), "Difference")  
df.comp$varNames<-as.factor(df.comp$varNames)

#ここまでで作図に必要なデータの加工が終了


####ここから図を作っていく

#theme set 真っ白なキャンバスを用意
ggplot()+theme_set(theme_bw(base_size=18,base_family="HiraKakuProN-W3"))

#plot
g <- ggplot(df.comp,aes(x = value,y = ..density..,colour=varNames,fill=varNames))
g <- g + geom_histogram(position = "identity",alpha = 0.5)
g <- g + geom_density(stat = "density",position = "identity",alpha = 0.4)
#facet
g<-g+facet_wrap(~parameter,scales="free")
g<-g + scale_color_brewer(palette="Set2")
g<-g + scale_fill_brewer(palette="Set2")
g<-g+ggtitle("Posterior Distributions")

#axis
g <- g + theme(axis.title.x = element_text(size=14), #axis title letter size
               axis.title.y = element_text(size=14)) 
g <- g + theme(axis.text.x = element_text(size=8),  #axis text letter size 
               axis.text.y = element_text(size=8))   
g <- g + theme(legend.title = element_text(size=14), #legend size
               legend.text = element_text(size=10),  #legend text size 
               legend.key= element_rect(colour = "white")) # legend mark-box colour



#MAP estimate-------------------------------------------------------
MAPs  = c(rep(NA,c((length(parameters2)-1)*(length(parameters2)-1)+1)))
for (i in 1:c((length(parameters2)-1)*(length(parameters2)-1)+1)){
    MAPs[i]    <- density(df.ext[,i])$x[which(density(df.ext[,i])$y == max(density(df.ext[,i])$y))]
}
names(MAPs) <- c(colnames(df.ext))


#情報量規準WAICとLOOの算出
#WAIC and Leave one out cross validation--------------------------------------------
WAIC<-c()
for (i in 1:Nvar){
    res<-waic(extract_log_lik(fit)[,c(Nsub*(i-1)+1):c(i*Nsub)])
    WAIC[i]<-res$waic
}

#Leave-One-Out cross validation
LOO<-c()
for (i in 1:Nvar){
    res2<-loo(extract_log_lik(fit)[,c(Nsub*(i-1)+1):c(i*Nsub)])
    LOO[i]<-res2$looic
}


#Final----------------------------------------------------------
# Print and Writing Results
AnalysisTime <- Sys.time() 
Results<-list(AnalysisTime=AnalysisTime,
              iterations=iterations,chains=chains,warmup=warmup,thin=thin,
              varNames=varNames,Summary=Summary,MAPs=MAPs)
InformationCriteria<-list(WAIC=WAIC,LOO=LOO)

print(Results)
print(InformationCriteria)
g

#writing 
ggsave("GaussiamDifference.png",g)

Written on January 11, 2017