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

# 多峰性とは

ヒストグラムを作ったときに、よく山が複数できることがあります。これを**多峰性(multimodal)**と呼びます。
以下に、３つの正規分布を結合させて３つの山を持つデータをつくってみます。

```R:
# 正規分布を３つ結合させる。
x <- c(rnorm(100, mean = 4, sd = 1),
       rnorm(100, mean = 8, sd = 1),
       rnorm(100, mean = 12, sd = 1))

# 階級幅1でヒストグラムを描く。
hist(x, breaks=seq(0,16,1))
```
![hist1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/424524/d2d793f7-bdd0-e24a-5766-ef218d4fd6af.png)


実際のデータでこの多峰性が現れた場合、データに複数の集団が含まれている可能性があります。しかし、ヒストグラムは階級幅の大きさによって形が変わってしまうため、本当の山(mode)の数はわかりません。例えば、先ほどのデータに対して階級幅を変えてヒストグラムを作ってみると、山の数が2個に見えるようになります。

```r:
# 階級幅2.1でヒストグラムを描く。
hist(x, breaks=seq(0,16,2.1))
```
![hist2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/424524/e44fa22d-346a-3365-205c-977ca655abf5.png)
この山の数を統計的な検定で決定したいと調べてみたところ、以下の資料に出会いました。

[ノンパラメトリックな多峰性検定-Silvermanの検定-とその古生物学への導入](https://www.jstage.jst.go.jp/article/kaseki/97/0/97_KJ00009874190/_article/-char/ja/)

資料では、多峰性検定の一つであるSilvermanの検定をとてもわかりやすく説明しており、とても勉強になりました。Rのパッケージも紹介していたのですが、見つけることができなかったので、以下のパッケージを試してみました。


[jenzoprさんのsilvermantestパッケージ（GitHub）](https://github.com/jenzopr/silvermantest)




# silvermantestパッケージの紹介

まず、導入してみます。

```r:
library(silvermantest)
```

`silverman.test`で検定を行うことができます。mode数が1のときは、キャリブレーションが必要となるので、`adjust=TRUE`とします。以下に、先ほどの３つの正規分布の結合データを検定してみました。

```r:

# mode数1で検定
silverman.test(x, 1, adjust = TRUE)

# mode数2で検定
silverman.test(x, 2)

# mode数3で検定
silverman.test(x, 3)
```

結果は以下のとおりです。mode数3でp値が0.05を超えました。mode数を3と考えて良さそうです。

```r:
Silvermantest: Testing the null hypothesis that the number of modes is <=  1 
The resulting p-value is  -4.5e-05 

Silvermantest: Testing the null hypothesis that the number of modes is <=  2 
The resulting p-value is  0.001001001 

Silvermantest: Testing the null hypothesis that the number of modes is <=  3 
The resulting p-value is  0.9479479 

Silvermantest: Testing the null hypothesis that the number of modes is <=  4 
The resulting p-value is  0.7767768 
```

上記の操作をまとめて行えるのが、`silverman.plot`です。mode数ごとに検定を行って、p値をグラフにしてくれます。グラフの赤い線が0.05のラインです。p値が0.05を初めて超えたmode数を適切と考えてよいようです。

```r:
silverman.plot(x)
```
![plot1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/424524/ac7886c8-fb5c-9ec5-7392-00b8e1425ec6.png)




さらに`densities.plot`で各mode数でのカーネル密度推定量のグラフを描くことができます。赤い線の部分が極大値を表しています。

```r:
densities.plot(x)
```

![plot2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/424524/ac59f0f8-548c-834f-73e2-3be7fd888530.png)


`densities.plot`のコードを加工して、ヒストグラムとカーネル密度推定量を重ねる関数を作ってみました。

```r:
# ヒストグラムとカーネル密度推定量を重ねる関数
hist_density = function(x, mode){
  
  # 極大値のインデックスを返す関数
  mode_ind = function(x, b) {
    y <- density(x, bw = b)$y
    d1 <- diff(y)
    signs <- diff(d1 / abs(d1))
    ind <- which(signs == -2)
    return(ind)
  }
  
  # 臨界バンド幅でのカーネル密度推定を行う。
  h0 <- h.crit(x, mode)
  d <- density(x, bw = h0)
  
  # グラフのx軸y軸の幅をを求める。
  y_max <- max(d$y)
  x_max <- max(abs(d$x))
  x_min <- min(abs(d$x))
  
  # カーネル密度推定量をプロット
  plot(d,
       main = paste("modes:", mode),
       xlim = c(x_min, x_max),
       ylim = c(0, y_max+0.1)
  )
  
  # 極大値のインデックスを求める。
  ind <- mode_ind(x, h0)
  
  x_values <- d$x
  y_values <- d$y
  
  # 極大値を赤くプロット。
  for (j in 1:mode){
    lines(c(x_values[ind[j]], x_values[ind[j]]),
          c(y_values[ind[j]] - 0.025, y_values[ind[j]] + 0.025),
          col = "red"
    )
  }
  
  # グラフを上書き。
  par(new=T)
  
  # ヒストグラムを描く
  hist(x,
       xlim = c(x_min, x_max),
       ylim = c(0, y_max+0.1),
       prob=T,
       ann=F)  

}
```

mode数3で実行してみます。

```r:
hist_density(x,3)
```
![hist3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/424524/a652c0eb-8647-e371-3c52-88688c1e1baf.png)


# Rの多峰性に関するパッケージ
　今回はSilvermanの検定とパッケージについて紹介しましたが、他にも多峰性に関する検定はたくさんあります。また、CRANには多峰性に関するパッケージである`multimode`があり、多くの検定だけでなく、mode tree や mode forest 、SiZer map などを扱うことができます。