非線形の相関の有無を検定できるHSIC testを試してみる

# はじめに
過去に書いた[この記事](https://qiita.com/purple_jp/items/69b97a76cfbb45048b63)で、いろいろな相関の指標を使ってみた経緯があります。
いろいろデータ分析をしていると、非線形な相関もうまく抽出したいというのがあります。MICもその指針となるのですが、他のアプローチもないものかと調べてみると、カーネル法を応用したHSIC testというものがあるらしいとのことで、これの使い方を調べて試した結果をメモしたくらいの備忘録です。

# HSIC testについて
再生核ヒルベルト空間を用いたノンパラメトリックな独立性検定のこと（Hilbert-Schmidt Independence Criterion test）。が、悲しいことに、私は、この辺の数理の理解はできてなく、時間を要しそうなので、参考になりそうなリンクだけ。
- [福水健次, 再生核ヒルベルト空間を用いた非線形データ解析法](https://www.ism.ac.jp/~fukumizu/papers/Nagoya_Nov2005.pdf)
- [福水健次, カーネル法入門 6.カーネル平均を用いたノンパラメトリック推論](https://www.ism.ac.jp/~fukumizu/OsakaU2014/OsakaU_6kernelMean.pdf)

# ひとまずRで試してみる
普通のHSICも使えるパッケージがありますが、純粋に拡張した[dHSIC(d variable HSIC)パッケージ](https://cran.r-project.org/web/packages/dHSIC/dHSIC.pdf))を使うのが良さそうです。HSICは2変数しか対応していませんでしたが、こちらは多変数が使える模様。

2変数で使う時は、dhsic.test(x, y)$p.valueでP値を求めることができる。帰無仮説は「2変数には(非)線形関係は無い」なので、小さければそれが棄却されるので、非線形関係も含めて相関があると言える。なお、[「Toru Imai, 21世紀の手法対決 (MIC vs HSIC)」](https://www.slideshare.net/motivic/tokyo-r-lt-25759212)によると、ノイズに強くても関係ありと判定しやすい模様。この点も試してみる。

## 線形
基本的には、ここに限らず、全て[過去記事](https://qiita.com/purple_jp/items/69b97a76cfbb45048b63)＋dhsic.testをやってみただけです。
あと、前よりノイズを増大させてみたものもあります。

![000011.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187459/a4e09f00-26ca-cc83-d613-f35f20a1ee90.png)

```r
library(dHSIC)
library(minerva)
set.seed(1)
x1 <- rnorm(n=100)
set.seed(2)
e <- rnorm(n=100)
y1 <- 4 * x1 + 5 + 2*e
plot(x1, y1)
cor(x1, y1, method="pearson")
cor(x1, y1, method="spearman")
mine(x1, y1)$MIC
dhsic.test(x1, y1)$p.value
```

```shell:出力結果
[1] 0.8159582 #pearson
[1] 0.8138854 #spearman
[1] 0.5862294 #MIC
[1] 0.000999001 #dHSIC p-value
```

## 3次曲線

![00000b.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187459/edb3c069-753d-f052-1b2c-8012b20da802.png)

```r
y2 <- 3*x1^3 + x1^2 + 2*x1+ 5 + 2*e
plot(x1, y2)
cor(x1, y2, method="pearson")
cor(x1, y2, method="spearman")
mine(x1, y2)$MIC
dhsic.test(x1, y2)$p.value
```

```shell:出力結果
[1] 0.8397053 #pearson
[1] 0.8214701 #spearman
[1] 0.6140474 #MIC
[1] 0.000999001 #dHSIC p-value
```

##2次曲線

![000007.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187459/316cffdb-4411-dc02-0f12-f0024a6b3daa.png)

```{r}
y3 <- 4*x1^2 + 5 + 3*e
plot(x1, y3)
cor(x1, y3, method="pearson")
cor(x1, y3, method="spearman")
mine(x1, y3)$MIC
dhsic.test(x1, y3)$p.value
```

```shell:出力結果
[1] 0.003019591 #pearson
[1] 0.01570957 #spearman
[1] 0.5286332 #MIC
[1] 0.000999001 #dHSIC p-value
```

## sinカーブ
dhsic.testのP値が高くなり帰無仮説を棄却できなくなりました。
周期性のあるデータの検出は苦手？この辺は要検討。

![000007.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187459/a2b6bf0f-4471-5532-b46f-2c0d6b34f8da.png)

```{r}
y4 <- 2*sin(5*x1) + e
plot(x1, y4)
cor(x1, y4, method="pearson")
cor(x1, y4, method="spearman")
mine(x1, y4)$MIC
dhsic.test(x1, y4)$p.value
```

```shell:出力結果
[1] -0.0958496 #pearson
[1] -0.08838884 #spearman
[1] 0.638089 #MIC
[1] 0.4785215 #dHSIC p-value
```

##乱数関係

![000005.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187459/3414ba5e-527a-54bd-1093-d66cb62b0c2c.png)

```{r}
set.seed(2)
y5 <- rnorm(n=100)
plot(x1, y5)
cor(x1, y5, method="pearson")
cor(x1, y5, method="spearman")
mine(x1, y5)$MIC
dhsic.test(x1, y5)$p.value
```

```shell:出力結果
[1] -0.1537208 #pearson
[1] -0.1537954 #spearman
[1] 0.2701778 #MIC
[1] 0.1488511 #dHSIC p-value
```

#最後に
sinカーブのdhsic.testで帰無仮説が棄却されなかった点は要考察。
数理面が分からないと、指標の限界も分からないので、要勉強ですね。。（当たり前）

あと、HSICは他の応用例も見かけるので、こんな感じで使ってみて終わりでなく、非線形の取り扱いの基礎のカーネル法やHSICはしっかりと理解して、適切な活用をしたいところです。

# 参考文献
- [R Package ‘dHSIC’](https://cran.r-project.org/web/packages/dHSIC/dHSIC.pdf)
- [福水健次, 再生核ヒルベルト空間を用いた非線形データ解析法](https://www.ism.ac.jp/~fukumizu/papers/Nagoya_Nov2005.pdf)
- [福水健次, カーネル法入門 6.カーネル平均を用いたノンパラメトリック推論](https://www.ism.ac.jp/~fukumizu/OsakaU2014/OsakaU_6kernelMean.pdf)
- [Toru Imai, 21世紀の手法対決 (MIC vs HSIC)](https://www.slideshare.net/motivic/tokyo-r-lt-25759212)

# 実行環境
- MacOS 10.15.5（19F101）
- R version 4.0.0 (2020-04-24) -- "Arbor Day"
- RStudio Version 1.3.959
- dHSIC 2.1
- minerva 1.5.8