t分布、Cauchy分布の区間と形状を視覚的につかむ

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

 


Figure. 左列: t分布、中央列: Cauchy分布、右列: 正規分布
それぞれ、自由度、位置パラメータ、標準偏差を1から5(上から下)に変化させている。

はじめに

上の図の左列のt分布。95%区間を示す縦の垂直線が斜めに見えません?これを「t錯視」といま適当に名付けました。

さて、今回は、 関東の有望なる “Mr.R & Stan” とでも言うべき彼に、t分布感と、それを感じるためのコード教わったので、少しやってみたいと思います。

Rでt分布とCauchy分布の乱数を出してみて、正規分布と比べてみます。

”いま、ここ”で、tとCauchyに対する新たな”気づき”を得ましょう。

ベイジアンにとっては、事前分布設定への新たな道が拓けるやもしれません。

本稿の目的

  1. t分布の乱数生成と可視化。自由度を変えてみる。
  2. Cauchy分布の乱数生成と可視化。尺度パラメータを変えてみる。
  3. 正規分布の乱数生成と可視化。分散パラメータ(標準偏差)を変えてみる。

下準備

デフォルトでは、t分布の乱数生成に際して平均と分散パラメータが指定できませんので、 {extraDistr},{metRorogy}パッケージを読み込みます。

###########################################################
#   Rデフォルトから、確率分布の増強と可視化 reported by D.H
###########################################################
invisible({rm(list=ls());gc();gc()})

#
#   目的: デフォルトにない特殊ケースでの乱数生成(t分布の平均と分散など利用)
#       1. t分布の自由度を変えて可視化する
#       2. Cauchy分布のScaleを変えて可視化する
#       3. 正規分布のSigmaを変えて可視化する
#

#初期設定-----------------------------------------------------
#分布を何個並べたいか(t分布は自由度、cauchyはScale、NormalはSigmaが変わる)
N=5
#何倍で変化させたいか(定数の設定)
Const=1
#乱数の数
R=1000
#-------------------------------------------------------------

#描画と加工のパッケージ-----------------------
library(ggplot2)
ggplot()+theme_set(theme_bw(base_size = 14,base_family = "HiraKakuProN-W3"))
library(Rmisc) #multiplot()
library(reshape2) #melt()
library(tidyverse)

#分布の追加パッケージ-----------------------------------
#install.packages("extraDistr")
library(extraDistr)
#install.packages("metRology")
library(metRology)

Studentのt分布”感”

まずはt分布を感覚的に把握していきます。

Stanで実装するとき、t分布には、(1)自由度nu、(2)平均(期待値)mu、(3)分散sigmaの3つのパラメータを使用します。 自由度はdfやk,あるいはnで表現されることもあります。本によっては、”正規化パラメータ”と呼ばれることもあります。

なお、本来的にt分布は自由度nuのみをパラメータとします。数学的な意味で。

平均はnuが1より大きいとき、分散はnuが2より大きい時のみ定義されます。

では、t分布の乱数を生成していきます。ここでは自由度を1から5まで変化させていきます。平均は0、標準偏差を5で固定します。

t分布はノイズに強い分布と言われますが、どのような特徴を持っているのでしょうか。

#Student_t Distribution------------------------------------------

#入れ物の用意
Student_t<-matrix(NA,N,R)
ci<-matrix(NA,N,2)
name_t<-c()

#自由度を変化させていく
for(i in 1: N){
    Student_t[i,]<-rt.scaled(R,  #1000個の乱数
             df=c(i*Const), #自由度は1から定数倍してN個
             mean=0, #平均0
            sd=5#標準偏差5
    )
    ci[i,]<-c(qt.scaled(c(0.025,0.975),#95%quantile
                      df=c(i*Const),
                      mean=0, 
                      sd=5)
              )#10 +- 1.96*sd
    name_t[i]<-sprintf("t:M=0,SD=5,DF=%.0f",i*Const)
    }
#乱数と95%HDIをtidyなデータフレームに
Student_t<-data.frame(Student_t)
Student_t$ID<-name_t
Student_t<-melt(Student_t,id=c("ID"))
CI<-data.frame(ci)
CI$ID<-name_t

#95%区間と乱数をIDで結合(merge)
df.t<-merge(Student_t,CI,by= c("ID"),all = T) %>%
    dplyr::select(ID,value,X1,X2)%>%
    dplyr::rename(Lower=X1,Upper=X2)

ここまでで、乱数生成と95%区間を得ることができました。 ggplotで描画します。

#plot
g1<-ggplot(df.t)+
    geom_density(aes(x=value,fill=ID),alpha=0.5)+
    geom_vline(aes(xintercept=Lower),linetype=3,size=0.8)+
    geom_vline(aes(xintercept=Upper),linetype=3,size=0.8)+
    scale_y_continuous(expand=c(0,0))+
    facet_wrap(~ID,nrow  = N)+
    scale_fill_brewer(palette="Greens")+
    xlim(-100,100)+ theme(legend.position="none")
plot(g1)

よいしょ。

自由度が1だと、とんでもない95%区間ですね。

自由度が大きくなるにつれて、裾が軽い(すぐゼロにへばりつく)分布になっています。

自由度が無限で正規分布に近似するそうですが、5でも結構正規分布に近い感じがします。

もっと、100000くらいじゃないと正規分布っぽくならないとおもってました。

Cauchy分布の分布感

Cauchy分布は、正規分布の分散パラメータ(標準偏差)の事前分布によく使用されます。

t分布の特殊系で、自由度1のt分布と同一です。

Cauchy分布のパラメータはLocation(位置)とScale(尺度)の二つ。平均は定義されていません。

階層線形モデルの傾きの階層パラメータなんかにも使われます(Gelmanによれば)。

弱情報事前分布(weakly informative prior)として、 ~ cauchy(0,2.5)という表記をよく見かけますが、あれってどんな形なんでしょうね。

位置パラメータを0、尺度パラメータを1から5まで変えて見てみましょう。

#Cauchy Distribution----------------------------------------
cauchy<-matrix(NA,N,R)
ci2<-matrix(NA,N,2)
name_cauchy<-c()
#Scaleを変化させていく
for(i in 1: N){
    cauchy[i,]<-rcauchy(R,  #1000個の乱数
                        location=0,#位置パラメータ
                        scale=i*Const#尺度パラメータ
    )
    ci2[i,]<-c(qcauchy(c(0.025,0.975),#95%quantile
                        location=0, 
                        scale=i*Const)
    )#10 +- 1.96*sd
    name_cauchy[i]<-sprintf("Cauchy:L=0,S=%.0f",i)
}
#乱数と95%HDIをtidyなデータフレームに
Cauchy<-data.frame(cauchy)
Cauchy$ID<-name_cauchy
Cauchy<-melt(Cauchy,id=c("ID"))
CI2<-data.frame(ci2)
CI2$ID<-name_cauchy

#95%区間と乱数をIDで結合(merge)
df.c<-merge(Cauchy,CI2,by= c("ID"),all = T) %>%
    dplyr::select(ID,value,X1,X2)%>%
    dplyr::rename(Lower=X1,Upper=X2)

#plot
g2<-ggplot(df.c)+
    geom_density(aes(x=value,fill=ID),alpha=0.5)+
    geom_vline(aes(xintercept=Lower),linetype=3,size=0.8)+
    geom_vline(aes(xintercept=Upper),linetype=3,size=0.8)+
    scale_y_continuous(expand=c(0,0))+
    facet_wrap(~ID,nrow  = N)+
    scale_fill_brewer(palette="Blues")+
    xlim(-100,100)+ theme(legend.position="none")
plot(g2)

せいやっ

うわもう外れ値みたいな値を出しまくり。

Cauchy分布のScaleパラメータというのは、正規分布の標準偏差みたいな感覚でいるとエライ目見ますね。

Scale=5で、その95%区間が±60を超えています。正規分布だとありえない。

でも、そのわりには、分布の盛り上がりは位置パラメータの0付近でしっかり出ています。

外れ値に強そう、というより、強すぎる感がありますね。

正規分布を比較として

正規分布についても、同じように出してみます。

あくまで比較のため、見にくいのですが、x軸はt分布およびCauchy分布に合わせます。

正規分布は平均0、標準偏差を1から5まで変えていきます。

#Normal Distribution----------------------------------------
Normal<-matrix(NA,N,R)
ci3<-matrix(NA,N,2)
name_Normal<-c()
#標準偏差を変化させていく
for(i in 1: N){
    Normal[i,]<-rnorm(R,0,i*Const)
    ci3[i,]<-c(qnorm(c(0.025,0.975),0,i*Const))
    name_Normal[i]<-sprintf("Norm:mu=0,sig=%.0f",i)
}
#乱数と95%HDIをtidyなデータフレームに
Normal<-data.frame(Normal)
Normal$ID<-name_Normal
Normal<-melt(Normal,id=c("ID"))
CI3<-data.frame(ci3)
CI3$ID<-name_Normal

#95%区間と乱数をIDで結合(merge)
df.n<-merge(Normal,CI3,by= c("ID"),all = T) %>%
    dplyr::select(ID,value,X1,X2)%>%
    dplyr::rename(Lower=X1,Upper=X2)

#plot
g3<-ggplot(df.n)+
    geom_density(aes(x=value,fill=ID),alpha=0.5)+
    geom_vline(aes(xintercept=Lower),linetype=3,size=0.8)+
    geom_vline(aes(xintercept=Upper),linetype=3,size=0.8)+
    scale_y_continuous(expand=c(0,0))+
    facet_wrap(~ID,nrow  = N)+
    scale_fill_brewer(palette="PuRd")+
    xlim(-100,100)+ theme(legend.position="none")+ theme(axis.title.y = element_text(size=0))
plot(g3)

おりゃ

ちっちゃ。区間せまっ。

と思われるかもしれませんが、Cauchyがひろすぎるのです。

3つを並べて見てみる。

このt-Cauchy-Normalが似た形に見えて実はぜんぜん違う”感”を味わいたいので、並べてみましょう。

(分布間の関係を上手く順序化してない。つまり対応がとれてないので、厳密な比較はしにくいと思います。すみません)

#MultiPlot
multiplot(g1,g2,g3,cols = 3)

<結論>

”事前分布はNormalでいいや”とか思うのはやめました。

データを見て、考えて、調べて、適切な事前分布を設定しようという姿勢が大事ですね。

Enjoy!!

Written on February 5, 2017