Follow @data_no_memo

メモ

個人的なメモです。他者にわかりやすく書くよりも未来の自分にわかりやすく書いています。なお、記事内容の正確さは保証できません。勉強中の身ですので、間違い等ご指摘頂けたら幸いです。

カーネル密度推定・Shinyを使ってみた

はじめに

 今回の記事は前半と後半に分かれる。前半でカーネル密度推定の概要を書く。後半でカーネル密度推定をコンテンツとしてShinyで簡単なweb アプリを作成する。はじめに両者の簡単な紹介をする。
 カーネル密度推定とは手元にあるデータからそのデータが従う確率分布を推定する方法である。これを用いれば、ヒストグラムから簡単に近似的な確率分布を求める事ができる。そうする事で視覚的に綺麗なデータを見せたりする事が可能となる。例えば、以下のようなヒストグラムがあったとしよう。

f:id:abcxyzonetwothree:20190601152039p:plain

 このヒストグラムからデータが得られた確率分布を表示したい時に、カーネル密度推定が使用される。それは以下のようになる。

f:id:abcxyzonetwothree:20190601152216p:plain

 しかしこの時、カーネル密度推定には分析者が決定するあるパラメータが存在し、そのパラメータの値によってはその近似された確率分布の見え方が異なってくる。例えば、そのパラメータを少しいじると以下のようになったりする。

f:id:abcxyzonetwothree:20190601152515p:plainf:id:abcxyzonetwothree:20190601152548p:plain

 このように、パラメータの値によっては近似される確率分布が粗くなったり、細かくなったりするのである。今回の記事では、このあたりも含めてカーネル密度推定では何が行われているのかを概説する。
 次に、Shinyについて簡単に説明する。ShinyはRのライブラリの1つであり、これを用いる事で簡単にwebアプリを作る事ができる。webアプリを使用すればデータの見せ方の選択肢が広がる。単に分析者が分析してその結果をポンと見せるのではなく、オーディエンスの入力に対してデータが動くというリアクティブなデータの見せ方が可能となる。例えば、オーディエンスが変数xと変数yをプロットしたいとすれば、任意の2変数を選択すればプロット図が出力されるようなwebアプリを作成しておけば、それを実現できる。
 今回の記事は以上2つを組み合わせる。すなわち、カーネル密度推定のパラメータをいじる事ができるwebアプリをShinyで作成する事が本記事の最終目標である。そのために、まずはカーネル密度推定についてもう少し詳しくみていく。

カーネル密度推定

 上記の確率分布はカーネル密度推定で求めた。それでは具体的に、どのような方法で密度を推定しているのか。なお、話を簡単にするために1次元の話に限定する。 ある確率分布p(x)からデータベクトルxが得られたとする。その中のデータのx'が領域Rに含まれる確率Pは以下のように表現できる。

$$P = \int_R p(x')dx' \tag{1}$$

 仮に、観測数Nがめちゃくちゃ大きければ(N→\inftyPは、Rに含まれるxの個数をkとすると、以下のように表現できる。

$$P = \frac{k}{N} \tag{2}$$

 今度は逆に、Rがめちゃくちゃ小さければ領域Rの体積をVとすると以下のように表現できる。

$$\int_R p(x')dx = p(x) V \tag{3}$$

(1)~(3)より以下のようになる。

$$p(x) = \frac{k}{NV} \tag{4}$$

 今、手元にあるデータの数は当然分かるのでNは既知である。一方で、kおよびVについてはわからない。そこで、それらを決める必要がある。
 まず、  Vについて。仮に、確率分布p(x)がある点を中心とする1辺の長さが hの立方体を領域RとすればV = h^{3}となる。一般化すると次元数 dの時、 V=h^{d}となる。今回は1次元を考えているので V = hとなる。すなわち、単なる線となる。この hは分析者が恣意的に決定する必要があるパラメータである。
 次に、kについて。Rとしてある点xを中心とする線を考えると、Rに入る標本数は以下のように表現できる。

$$k = \sum_{i+1}^N W(\frac{x - x_i}{h})\tag{5}$$

ただし、

$$W(x)=1   (if   |x| < \frac{1}{2}),  0 (else)$$

である。
 本当にこの W(x)によって Rに入る標本数を数え挙げられているのか。ここで、具体的に考えてみる。ここで x = (1, 2, 2, 3, 4, 4, 4)を考える。h  = 3とし1を中心とすると、1から +1.5 -1.5の範囲に収まる値の数だからk= 3となる。これを上式で当てはめて考えると以下のようになる。

$$k =W(\frac{1 - 1}{3}) + W(\frac{1 - 2}{3}) + W(\frac{1 - 2}{3}) + W(\frac{1 - 3}{3}) + W(\frac{1 - 4}{3}) + W(\frac{1 - 4}{3}) + W(\frac{1 - 4}{3})$$ $$=W(0) + W(-\frac{1}{3}) + W(-\frac{1}{3}) + W(-\frac{2}{3}) + W(-1) + W(-1) + W(-1)$$ $$=1 + 1 + 1 + 0 + 0 + 0 = 3$$

このように、ちゃんとk = 3となった。 ここで、(5)式と V = hを(4)式に代入すると以下のようになる。

$$p(x) = \frac{1}{Nh} \sum_{i+1}^N W(\frac{x - x_i}{h})$$

 実は、関数W(x)は他にも色なものを当てはめる事ができる。ここでは非連続的な変数を考えたが標準正規分布で代用すれば連続的な確率分布に近似できる。実は、冒頭のカーネル密度推定は標準正規分布を用いて密度推定を行っている。ここでどのような W(x)を考えるかは分析者が恣意的に決定する必要がある。
 以上からカーネル密度推定ではh W(x)を分析者が決定する必要がある。しかし、これらを決めてやれば近似的な確率分布を図示化する事ができる。

Shinyでカーネル密度推定の様子

 それではカーネル密度推定の様子をShinyを用いてwebアプリにしてみる。先に言及したようにカーネル密度推定を行う際には2つの事を決定する必要がある。1つ目がhの値である。この値によって推定される確率分布の形が異なる。先に答え合わせをするとhの値が大きくなればなるほど推定される確率分布は粗くなり、小さくなればなるほど推定される確率分布は細かくなる。次に決めるものは W(x)である。これには標準正規分布などいくつかの確率分布が想定する事ができる。
 それでは、以上2つの事柄を自由に決定しそれに合わせて反応するwebアプリケーションを作る。なお、Shinyの扱い方については以下の記事などを参照されたい。

qiita.com

 まず元になるデータを発生させる。

set.seed(1234)
x1 <- rnorm(1000, mean = 50, sd = 10)
x2 <- rnorm(1000, mean = 10, sd = 20)
x <- c(x1, x2)

 次に必要なlibraryを読み込む。

library(shiny)

 次にwebアプリのui部分を作成する。上述したように、2つの事柄(hの値とW(x)の形)を決定できるようにしている。

ui <- fluidPage(
    
    titlePanel("カーネル密度推定"),
    
    sidebarLayout(
        sidebarPanel(
            sliderInput("slider", label = h3("hの値"),
                        min = 0.001, max =3, value = 1, step = 0.001),
            selectInput("Kernel", label = h3("W(x)の形"),
                        choices = list("gaussian" = "gaussian", 
                                       "epanechnikov" = "epanechnikov", 
                                       "rectangular" = "rectangular",
                                       "triangular" = "triangular",
                                       "biweight" = "biweight",
                                       "cosine" = "cosine",
                                       "optcosine" = "optcosine"),
                        selected = "gaussian")
        ),
        mainPanel(
            h3("近似された確率分布"),
            plotOutput("ggplot"),
            h3("元のヒストグラム"),
            plotOutput("gghistogram")
        )
    )
    
)

 次にserverの部分を以下のように作る。

server <- function(input, output) {
    output$ggplot <- renderPlot({
        library(tidyverse)
        i <- input$slider
        kernel <- input$Kernel
        ggplot(mapping = aes(x = x)) + 
            geom_density(adjust = i, kernel = kernel) + 
            theme_bw()
    })
    output$gghistogram <- renderPlot({
        library(tidyverse)
        i <- input$slider
        kernel <- input$Kernel
        ggplot(mapping = aes(x = x)) + 
            geom_histogram(binwidth = 2) + 
            theme_bw()
    })
}

 最後に結果を見る。

shinyApp(ui = ui, server = server)

 できたwebアプリケーションは以下のリンクから確認する事ができる。実際に試されたい。

date-no-memo.shinyapps.io

 念の為、そのwebアプリケーションの見た目は以下の通りである。

f:id:abcxyzonetwothree:20190602200006p:plain