カーネル密度推定・Shinyを使ってみた
はじめに
今回の記事は前半と後半に分かれる。前半でカーネル密度推定の概要を書く。後半でカーネル密度推定をコンテンツとしてShinyで簡単なweb アプリを作成する。はじめに両者の簡単な紹介をする。
カーネル密度推定とは手元にあるデータからそのデータが従う確率分布を推定する方法である。これを用いれば、ヒストグラムから簡単に近似的な確率分布を求める事ができる。そうする事で視覚的に綺麗なデータを見せたりする事が可能となる。例えば、以下のようなヒストグラムがあったとしよう。
このヒストグラムからデータが得られた確率分布を表示したい時に、カーネル密度推定が使用される。それは以下のようになる。
しかしこの時、カーネル密度推定には分析者が決定するあるパラメータが存在し、そのパラメータの値によってはその近似された確率分布の見え方が異なってくる。例えば、そのパラメータを少しいじると以下のようになったりする。
このように、パラメータの値によっては近似される確率分布が粗くなったり、細かくなったりするのである。今回の記事では、このあたりも含めてカーネル密度推定では何が行われているのかを概説する。
次に、Shinyについて簡単に説明する。ShinyはRのライブラリの1つであり、これを用いる事で簡単にwebアプリを作る事ができる。webアプリを使用すればデータの見せ方の選択肢が広がる。単に分析者が分析してその結果をポンと見せるのではなく、オーディエンスの入力に対してデータが動くというリアクティブなデータの見せ方が可能となる。例えば、オーディエンスが変数xと変数yをプロットしたいとすれば、任意の2変数を選択すればプロット図が出力されるようなwebアプリを作成しておけば、それを実現できる。
今回の記事は以上2つを組み合わせる。すなわち、カーネル密度推定のパラメータをいじる事ができるwebアプリをShinyで作成する事が本記事の最終目標である。そのために、まずはカーネル密度推定についてもう少し詳しくみていく。
カーネル密度推定
上記の確率分布はカーネル密度推定で求めた。それでは具体的に、どのような方法で密度を推定しているのか。なお、話を簡単にするために1次元の話に限定する。 ある確率分布からデータベクトルが得られたとする。その中のデータのが領域に含まれる確率は以下のように表現できる。
$$P = \int_R p(x')dx' \tag{1}$$
仮に、観測数がめちゃくちゃ大きければ()は、に含まれるの個数をとすると、以下のように表現できる。
$$P = \frac{k}{N} \tag{2}$$
今度は逆に、がめちゃくちゃ小さければ領域の体積をとすると以下のように表現できる。
$$\int_R p(x')dx = p(x) V \tag{3}$$
(1)~(3)より以下のようになる。
$$p(x) = \frac{k}{NV} \tag{4}$$
今、手元にあるデータの数は当然分かるのでは既知である。一方で、およびについてはわからない。そこで、それらを決める必要がある。
まず、について。仮に、確率分布がある点を中心とする1辺の長さがの立方体を領域とすればとなる。一般化すると次元数の時、となる。今回は1次元を考えているのでとなる。すなわち、単なる線となる。このは分析者が恣意的に決定する必要があるパラメータである。
次に、について。としてある点を中心とする線を考えると、に入る標本数は以下のように表現できる。
$$k = \sum_{i+1}^N W(\frac{x - x_i}{h})\tag{5}$$
ただし、
$$W(x)=1 (if |x| < \frac{1}{2}), 0 (else)$$
である。
本当にこのによってに入る標本数を数え挙げられているのか。ここで、具体的に考えてみる。ここでを考える。とし1を中心とすると、1から、の範囲に収まる値の数だからとなる。これを上式で当てはめて考えると以下のようになる。
$$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$$
このように、ちゃんととなった。 ここで、(5)式とを(4)式に代入すると以下のようになる。
$$p(x) = \frac{1}{Nh} \sum_{i+1}^N W(\frac{x - x_i}{h})$$
実は、関数は他にも色なものを当てはめる事ができる。ここでは非連続的な変数を考えたが標準正規分布で代用すれば連続的な確率分布に近似できる。実は、冒頭のカーネル密度推定は標準正規分布を用いて密度推定を行っている。ここでどのようなを考えるかは分析者が恣意的に決定する必要がある。
以上からカーネル密度推定ではとを分析者が決定する必要がある。しかし、これらを決めてやれば近似的な確率分布を図示化する事ができる。
Shinyでカーネル密度推定の様子
それではカーネル密度推定の様子をShinyを用いてwebアプリにしてみる。先に言及したようにカーネル密度推定を行う際には2つの事を決定する必要がある。1つ目がの値である。この値によって推定される確率分布の形が異なる。先に答え合わせをするとの値が大きくなればなるほど推定される確率分布は粗くなり、小さくなればなるほど推定される確率分布は細かくなる。次に決めるものはである。これには標準正規分布などいくつかの確率分布が想定する事ができる。
それでは、以上2つの事柄を自由に決定しそれに合わせて反応するwebアプリケーションを作る。なお、Shinyの扱い方については以下の記事などを参照されたい。
まず元になるデータを発生させる。
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つの事柄(の値と)を決定できるようにしている。
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アプリケーションは以下のリンクから確認する事ができる。実際に試されたい。
念の為、そのwebアプリケーションの見た目は以下の通りである。