multimodeパッケージで多峰性を視覚的に検討する

はじめに

　前回、多峰性検定であるSilvermanの検定について以下の記事を書きました。

　多峰性をSilvermanの検定で確かめる

　しかし、多峰性検定がうまく使えないことがあり、調べたところ以下の資料に出会いました。

　ノンパラメトリック・スムージング理論とその応用（樋田勉先生）

　この資料には、mode tree、mode forestなど多峰性を視覚的に検討する方法について書かれており、非常に勉強になりました。また、これらの方法を利用するにあたり、CRANにあるmultimodeパッケージがとても優れていましたので、紹介したいと思います。

データの準備

　以下のように、正規分布を３つ結合したデータを用意しました。ヒストグラムをみると、山(mode)を３つ確認できます。以降、このデータを使っていきます。

# ３つの正規分布を結合する。
x <- c(rnorm(100, mean = 5, sd = 1),
       rnorm(100, mean = 20, sd = 1),
       rnorm(100, mean = 100, sd = 1))

# ヒストグラムを作成。
hist(x, bins=1, breaks = seq(1,120,1))

multimodeの導入

multimodeは、CRANにありますので、以下のようにインストールできます。

install.packages("multimode")

他のパッケージと同様に、読み込みます。

library(multimode)

多峰性検定

　multimodeパッケージでは、modetestで複数の多峰性検定が利用できます。今回は、Silvermanの検定を使いたいと思います。以下のように、mode数ごとにp値を計算し、グラフにしました。

# mode数ごとにp値を計算する。
p_value = c()
for (k in 1:10){
  result <- modetest(data = x,
            mod0 = k,
            method = "SI",
            submethod = 2,
            B = 100)

  p_value <- append(p_value, result$p.value,)
}

# グラフを作成する。
plot(p_value,
     type = "o",
     xlab = "number of mode",
     ylab = "p value")

　上記のグラフを見てみると、p値が0.05を超えることがなく、mode数を確定できません。mode数を10まで増やして検定しても同様でした。多峰性検定では「分布が連続である」という仮定が必要で、今回のデータのように分布が離れている場合には適切ではなく、mode数を確定できないことが多いようです。
　このような場合に、視覚的に多峰性を検出する方法がmode tree, mode forest, SiZer map です。

mode tree

カーネル密度推定では、バンド幅を変えるとグラフの形が変わり、modeの位置も変化します。 このバンド幅とmodeの位置の変化を視覚化したものがmode treeです。multimodeパッケージでは、modetreeでmode treeを作成できます。以下に作ってみます。

# mode treeを作成する。
modetree(x)

　できたmode treeを見てみると、まず、13付近のmodeから100付近のmodeが分かれます。さらにバンド幅が小さくなると、13付近のmodeから5付近と20付近のmodeへ分かれています。このことからmode数は 2（13付近と100付近）、または 3（5付近と20付近と100付近）であることが示唆されます。mode数が 2 なのか 3 なのか、mode forestを使うと細かく検討ができます。

mode forest

mode forestは、標本からリサンプリングを行い、それぞれのリサンプリング標本からできたmode treeを重ねてつくられたものです。着色が濃いほど、modeの存在の可能性が高くなります。multimodeパッケージでは、modeforestでmode forestを作成できます。以下に作ってみます。

# mode forestを作成する。
modeforest(x)

　5付近、20付近、100付近が濃く、modeが存在する可能性が高いことを示しています。一方、13付近は薄く、modeが存在する可能性が低いことを示しています。このことからmode数は 3（5付近と20付近と100付近）であることが示唆されます。

SiZer map

　SiZer map は、mode treeやmode forestと同様にバンド幅を変えてカーネル密度推定を行っています。そして、各バンド幅でカーネル密度推定量が有意に増加しているところを青、減少しているところを赤、微分係数がゼロであるところを紫で表しています。したがって、青→(紫)→赤の順で変化している場所が、カーネル密度推定量の山のピーク(mode)を表しています。multimodeパッケージでは、sizerでSiZer mapを作成できます。以下に作ってみます。

# SiZer mapを作成する。
sizer(x)

5付近、20付近、100付近に青→紫→赤の変化が見られますので、modeもこの位置であると考えられます。mode treeを重ね描きすると、わかりやすくなります。以下に作ってみます。

# mode treeを重ね描きする。
modetree(x,
         logbw = TRUE,
         addplot = TRUE,
         col.lines = "white")

　青→紫→赤の変化の部分にmode treeが描かれており、modeの存在が示唆されています。

ノンパラメトリックとパラメトリックの多峰性

　今回の記事は、ノンパラメトリックな多峰性を扱いましたが、パラメトリックな多峰性についてはガウス混合分布を利用することができます。