目次

・母集団の分布がわからないときに分布を推定する

・カーネル密度推定

　　・シルバーマンの経験則

・Rで行うカーネル密度推定

　　・カーネル関数の影響

　　・母集団が正規分布以外のとき（ガンマ分布）

　　・母集団が正規分布以外のとき（二項分布）

・カーネル密度推定の応用：ヴァイオリンプロット

母集団の分布がわからないときに分布を推定する

　今まで、私たちはあらかじめ分布がわかっているデータについて解析を行ってきた。あるいは、こういうデータならこういう分布になるだろうと、ある程度のあたりをつけて解析を先に進めることもあるだろう。例えば、正の値しかとらない連続値データなら、ガンマ分布に従いそうだから、shapeパラメータやscaleパラメータをデータから推定し、当てはめるだろう。正の値しかとらない連続値データでも、左右対称な分布なようなら正規分布に従うと考えても差し支えないから、平均や分散パラメータを推定するだろう。しかしながら、実際にはデータから元の分布を推定することは困難が伴う。そこで、なるべく少ない仮定でデータが従いそうな分布を構築することを目指す。その時に用いられる手法が、カーネル密度推定kernel density estimationである。

カーネル密度推定

　カーネル密度推定は、サンプルサイズnの標本から、以下の式を用いて確率密度分布を推定するノンパラメトリックな手法である。

hはバンド幅と呼ばれ、その決め方は後述する。K(x)はカーネル関数と呼ばれ、標準ガウス関数とすることが多い。カーネル関数の効果も後程考えるとして、標準ガウス関数をカーネル関数とするときに、この式がどういう意味を持つかを考えよう。

　以下のように標準ガウス関数を代入する。

すると、この式は、xiを平均、hを標準偏差とする正規分布の確率密度関数を足し上げて、サンプルサイズnで割ったもの（平均）であると理解できる。一つ一つのデータにカーネル関数である正規分布を当てはめて、その和でデータの従う確率密度関数を表現しようというわけだ。ここで気になるのは、hの決め方である。hが大きければ全体的になだらかで連続的な分布になるけど、外れ値も出やすい分布になる。hか小さければ外れ値は出にくいけど、とびとびの分布になりやすい。

そういうわけで、適度なバンド幅hをこちらで決定する必要がある。

・シルバーマンの経験則

　バンド幅の決め方に、絶対というものはないが、標準ガウス関数を使うときには、しばしばシルバーマンの経験則Silverman's ‘rule of thumb’ (Silverman, 1986, Density Estimation) が用いられる。シルバーマンの経験則は、データの標準偏差、四分位数範囲、そしてサンプルサイズを用いてhを以下のように決める。

なあにこれぇ、って感じだけど、とりあえずこの式にのっとってhを決めると、経験則上は良い感じ！という式である。

　式としてはシンプルで、標準偏差か、四分位範囲/1.34のうち小さい方に0.9とサンプルサイズの-1/5乗をかける。正規分布において、平均を中心として50%の領域は、母標準偏差を約0.67倍した範囲であるから、四分位範囲の示す値はおおむね0.67×2 = 1.34σとなる。

ゆえに、データが正規分布に従い、サンプルサイズ十分大きければ、標準偏差と四分位範囲/1.34はともに似たような値になると期待される。サンプルサイズが小さいときほど、1つのデータの影響が大きく標準偏差は変動が大きいが、四分位範囲への影響は弱く、四分位範囲のほうが小さい値を示すことが多い。

　サンプルサイズが大きいとき、n^(-1/5)はより小さな値を示すようになる。ゆえに、サンプルサイズが大きいときほど、バンド幅は狭くなる。これは、サンプルサイズが大きいほど、サンプルの信頼度が上がり、より狭い正規分布に当てはめることにする、と考えることができるだろう。

Rで行うカーネル密度推定

　では、Rでカーネル密度推定を行ってみる。density関数を用いることで、簡単に行うことができる。平均10、標準偏差5の正規分布に従うデータを2個だけ生成してみよう。併せていくつかの作図を行う。

------------------------------------------------------

library(plotn)

n <- 2

m <- 10

s <- 5

d <- rnorm(n, mean = m, sd = s)#正規分布に従うデータを2個生成

k <- density(d)#カーネル密度推定

plotn(0, 0, xlim = range(k$x), 

      ylim = range(c(k$y, dnorm(k$x, mean = m, sd = s))), 

      xlab = "X", ylab = "density",

      col.dot = NA)#図1の描画

for(i in 1:n){

  overdraw(segments(d[i], 0, 

                    d[i], 0.05*max(c(k$y, dnorm(k$x, mean = m, sd = s))), 

                    col = "#000000"),#データ点の描画 

           points(k$x, (1/n)*dnorm(k$x, mean = d[i], sd = k$bw), 

　　　　　　　　　　　type = "l", col = "blue"))#データ点周りの正規分布の1/n

}

overdraw(points(k$x, k$y, type = "l"),#カーネル密度推定値

         points(k$x, dnorm(k$x, mean = m, sd = s), type = "l", col = "red"))#母集団

------------------------------------------------------

kはlist形式で、いくつかのデータが入っているが、重要なのはk$x、k$y、k$bwである。k$yはk$xに対するカーネル密度推定値である。k$bwはバンド幅である。図1のように、カーネル密度推定値は青で描かれたデータ点周りの正規分布の確率密度関数の和になっていることがわかるだろう。データがあまりに少ないので、母集団（赤線）とは似つかない形になっているが、データ点が増えれば問題がなくなることをここから確かめていく。

　ちなみに、デフォルトのバンド幅は以下のように求められる。確かに一致している。

------------------------------------------------------

k$bw#バンド幅

## [1] 0.9680107

sd(d)#標準偏差

## [1] 2.341334

(quantile(d)[4] - quantile(d)[2])/1.34#四分位範囲/1.34

##      75% 

## 1.235503

0.9 * min(sd(d), (quantile(d)[4] - quantile(d)[2])/1.34) * (n)^(-1/5)#バンド幅の手計算

## [1] 0.9680107

------------------------------------------------------

今度は、サンプルサイズを5まで増やしてカーネル密度推定値の出来上がっている様子を確認してみよう。作図も少し工夫を凝らす。

------------------------------------------------------

n <- 5

m <- 10

s <- 5

d <- rnorm(n, mean = m, sd = s)

d <- d[order(d)]

k <- density(d)

plotn(0, 0, xlim = range(k$x), 

      ylim = range(c(k$y, dnorm(k$x, mean = m, sd = s))), 

      xlab = "X", ylab = "density",

      col.dot = NA)#図2の描画

ks <- NULL

for(i in 1:n){

  ks <- cbind(ks, (1/n)*dnorm(k$x, mean = d[i], sd = k$bw))#各正規分布の確率密度を計算

}

for(i in 1:(n-1)){

  overdraw(polygon(c(k$x[1], k$x, k$x[length(k$x)]), 

                   c(0,rowSums(ks[,1:(n-i+1)]),0), 

                   col = .default_col[i], border = NA))#各正規分布の確率密度の和の描画

}

overdraw(polygon(c(k$x[1], k$x, k$x[length(k$x)]), 

                 c(0, ks[,1], 0), 

                 col = .default_col[n], border = NA))#最後の1つの正規分布の確率密度の和の描画

for(i in 1:n){

  overdraw(segments(d[i], 0, 

                    d[i], 0.05*max(c(k$y, dnorm(k$x, mean = m, sd = s))), 

                    col = "#000000"),#データ点の描画

           points(k$x, (1/n)*dnorm(k$x, mean = d[i], sd = k$bw), type = "l"))#データ点周りの正規分布の1/n

}

overdraw(points(k$x, k$y, type = "l", lwd = 2),#カーネル密度推定値

         points(k$x, dnorm(k$x, mean = m, sd = s), type = "l", col = "red"))#母集団

------------------------------------------------------

平均値の小さい方の正規分布から順に積み上げて重ねたような図にしてみた。カーネル密度推定値が、正規分布の和で出来上がっていることがよくわかるだろう。まだまだ、母集団には程遠いが、徐々に一山形になりつつあることがわかる。

　では、一気にサンプルサイズが10000まで行くとどうなるだろう。

------------------------------------------------------

n <- 10000

m <- 10

s <- 5

d <- rnorm(n, mean = m, sd = s)

k <- density(d)

plotn(0, 0, xlim = range(k$x), 

      ylim = range(c(k$y, dnorm(k$x, mean = m, sd = s))), 

      xlab = "X", ylab = "density",

      col.dot = NA)#図3の描画

for(i in 1:n){

  overdraw(segments(d[i], 0, 

                    d[i], 0.05*max(c(k$y, dnorm(k$x, mean = m, sd = s))), 

                    col = "#000000"))

}

overdraw(points(k$x, k$y, type = "l"),

         points(k$x, dnorm(k$x, mean = m, sd = s), type = "l", col = "red"))

------------------------------------------------------

すると、カーネル密度推定値は、母集団によく一致するようになることがわかる。

・カーネル関数の影響

　実は、デフォルトのカーネル関数は標準ガウス関数（Gaussian）であるが、他のカーネル関数を使うこともできる。例えば、一様関数や三角型関数である。

標準ガウス関数を用いた場合、すべての実数で定義されているから、あるデータ点は、遠くのデータ点の確率密度関数の影響をわずかながらに受ける。しかし、三角型関数や一様関数をカーネル関数にした場合、あるところ境界に遠くのデータ点の影響を受けなくなる。より「境界が明瞭」な分布を描きたいときはこれらのカーネル関数も選択肢にあがるが、そのかわり滑らかではない点も生ずる。サンプルサイズが十分大きければ、どのカーネル関数でもあまり変わらなくなる。

　三角型関数や一様関数を使ってカーネル密度推定を行うと以下の通り。

------------------------------------------------------

k <- density(d, kernel = "triangular")#カーネル関数を三角型関数に変更

plotn(0, 0, xlim = range(k$x), 

      ylim = range(c(k$y, dnorm(k$x, mean = m, sd = s))), 

      xlab = "X", ylab = "density",

      col.dot = NA)#図4の描画

for(i in 1:n){

  overdraw(segments(d[i], 0, 

                    d[i], 0.05*max(c(k$y, dnorm(k$x, mean = m, sd = s))), 

                    col = "#000000"))

}

overdraw(points(k$x, k$y, type = "l"),

         points(k$x, dnorm(k$x, mean = m, sd = s), type = "l", col = "red"))

k <- density(d, kernel = "rectangular")#カーネル関数を一様関数に変更

plotn(0, 0, xlim = range(k$x), 

      ylim = range(c(k$y, dnorm(k$x, mean = m, sd = s))), 

      xlab = "X", ylab = "density",

      col.dot = NA)#図5の描画

for(i in 1:n){

  overdraw(segments(d[i], 0, 

                    d[i], 0.05*max(c(k$y, dnorm(k$x, mean = m, sd = s))), 

                    col = "#000000"))

}

overdraw(points(k$x, k$y, type = "l"),

         points(k$x, dnorm(k$x, mean = m, sd = s), type = "l", col = "red"))

------------------------------------------------------

・母集団が正規分布以外のとき（ガンマ分布）

　母集団が正規分布なら、標準ガウス関数の和で母集団を表現できることは想像に難くない。では、母集団が正規分布ではなかったときも、カーネル密度推定はうまくいくだろうか。例えば、ガンマ分布は左右非対称で、0以下の値をとらない。以下で試してみよう。

------------------------------------------------------

n <- 10000

sh <- 3

sc <- 10

d <- rgamma(n, shape = sh, scale = sc)

k <- density(d)

plotn(0, 0, xlim = range(k$x), 

      ylim = range(c(k$y, dgamma(k$x, shape = sh, scale = sc))), 

      xlab = "X", ylab = "density",

      col.dot = NA)#図6の描画

for(i in 1:n){

  overdraw(segments(d[i], 0, 

                    d[i], 0.05*max(c(k$y, dgamma(k$x, shape = sh, scale = sc))), 

                    col = "#000000"))

}

overdraw(points(k$x, k$y, type = "l"),

         points(k$x, dgamma(k$x, shape = sh, scale = sc), type = "l", col = "red"))

------------------------------------------------------

すると、左右対称ではない様子はうまく表現できているように見える。標準ガウス関数をカーネル関数にする以上、どうしても0以下の確率密度が推定されるが、それでも全体的にうまく推定ができている。

・母集団が正規分布以外のとき（二項分布）

　離散分布にも適用できるかも気になるところだ。そこで、二項分布から生成されたデータをカーネル密度推定してみる。

------------------------------------------------------

n <- 10000

si <- 100

pr <- 19/20

d <- rbinom(n, size = si, prob = pr)

k <- density(d)

plotn(0, 0, xlim = range(k$x), 

      ylim = range(c(k$y, dbinom(0:si, size = si, prob = pr))), 

      xlab = "X", ylab = "density",

      col.dot = NA)#図7の描画

for(i in 1:n){

  overdraw(segments(d[i], 0, 

                    d[i], 0.05*max(c(k$y, dbinom(0:si, size = si, prob = pr))), 

                    col = "#00000003",

                    lwd = 10))

}

overdraw(points(k$x, k$y, type = "l"),

         points(0:100, dbinom(0:si, size = si, prob = pr), type = "l", col = "red"),

         points(0:100, dbinom(0:si, size = si, prob = pr), col = "red"))

------------------------------------------------------

まあ、そもそも離散分布に正規分布の当てはめを行おうという、無理やりなことをしているわけだから歪な分布になる。本来は存在しない、非整数値にも確率密度が推定されることになる。おそらくは、各整数xがとる確率をx - 1/2からx + 1/2の範囲の面積とすれば、二項分布の確率質量と近い値になるだろう。上記のように、正規近似が成り立たなそうな状況では歪な分布になるが、より左右対称な状況では、もっとましになる。

------------------------------------------------------

n <- 10000

si <- 100

pr <- 5/6

d <- rbinom(n, size = si, prob = pr)

k <- density(d)

plotn(0, 0, xlim = range(k$x), 

      ylim = range(c(k$y, dbinom(0:si, size = si, prob = pr))), 

      xlab = "X", ylab = "density",

      col.dot = NA)#図8の描画

for(i in 1:n){

  overdraw(segments(d[i], 0, 

                    d[i], 0.05*max(c(k$y, dbinom(0:si, size = si, prob = pr))), 

                    col = "#00000003",

                    lwd = 10))

}

overdraw(points(k$x, k$y, type = "l"),

         points(0:100, dbinom(0:si, size = si, prob = pr), type = "l", col = "red"),

         points(0:100, dbinom(0:si, size = si, prob = pr), col = "red"))

------------------------------------------------------

カーネル密度推定の応用：ヴァイオリンプロット

　カーネル密度推定の応用として、ヴァイオリンプロットviolin plotと呼ばれる作図がある。例えば、箱ひげ図とともに、カーネル密度推定を作図する（あるいは箱ひげを書かない場合もある）。その形がヴァイオリンみたいなので、こう呼ばれる。

　ヴァイオリンプロットを行うパッケージはいくつかあるが、私が作ったplotnパッケージでもvioplotn関数によって作図できる。例えば、以下は一般的な箱ひげ図である。

------------------------------------------------------

n <- 20

y1 <- rnorm(n, mean = 10, sd = 5)

y2 <- rnorm(n, mean = 15, sd = 2)

d <- data.frame(group = rep(c("A","B"), each = n), y = c(y1, y2))

boxplotn(y ~ group, data = d, jitter.method = "swarm")#図9の描画

------------------------------------------------------

対して、ヴァイオリンプロットは以下のような感じだ。

------------------------------------------------------

vioplotn(y ~ group, data = d, jitter.method = "swarm", trim = T)#図10の描画

------------------------------------------------------

このような作図によって、大まかなデータの分布をさらに視覚的に示すことが可能となる。ちなみに、上記では、trim引数をTとすることで、実際に得られたデータよりも大きい値や小さい値の密度推定を描かないようにしている。実際に得られたデータ外は、外挿となるため、各日にわかりそうなことだけを描く措置である。trim = Fとすれば、滑らかな曲線がデータが描かれる。