非線形な関係も扱える相関係数MICとヒルベルト・シュミット独立性検定HSIC

変数の独立性を測る指標として有名なものにピアソンの相関係数(Pearson's correlation coefficient)があります。ですがそれはあくまで線形な（直線的な）関係を扱う指標で、たとえば二次関数のようなU字型の関係があるときには「相関がない」と判定してしまいます。このような非線形な関係を取り扱う指標として、MICとHSICが知られているので、そのPythonパッケージを使ってみました。

以下のコードは全て Google Colaboratory で動作確認済みです。

# MIC (Maximal Information Coefficient) を計算するPythonパッケージ
次のようにして pip install できます。

```Python
!pip install minepy
```

インストールできたら、次のようにして準備しておきます。

```Python
from minepy import MINE
mine = MINE()
```

# HSIC (Hilbert-Schmidt Independence Criterion) を計算するPythonパッケージ
次のように git clone します。

```Python
!git clone https://github.com/amber0309/HSIC.git
```

次のようにして準備しておきます。

```Python
from HSIC.HSIC import hsic_gam
```

# 一次関数
まずは一次関数について調べてみます。

```Python
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

X = np.linspace(-1, 1, 51)
Y = X
plt.scatter(X, Y)
```
![out6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/39d8c80b-b13e-4fb3-3b27-e1627fea2b3f.png)

ピアソンの相関係数は

```Python
np.corrcoef(X, Y)[0, 1]
```
```
1.0
```

MICは

```Python
mine.compute_score(X, Y)
mine.mic()
```

```
0.9999999999999998
```

HSSCは、次の値が大きいほど独立性が低く、小さいほど独立性が高いそうです。

```Python
testStat, thresh = hsic_gam(X.reshape(len(X), 1), Y.reshape(len(Y), 1))
testStat / thresh
```

```
22.610487654066624
```

一次関数にノイズを加えてみましょう。

```Python
import numpy as np

X = np.linspace(-1, 1, 51)
Y = X + np.random.rand(51,)
plt.scatter(X, Y)
```
![out10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/fb0f90d0-1cc8-9322-710c-b324a8f600d2.png)

ピアソンの相関係数は

```Python
np.corrcoef(X, Y)[0, 1]
```

```
0.9018483832673769
```

MICは

```Python
mine.compute_score(X, Y)
mine.mic()
```

```
0.9997226475394071
```

HSICは

```Python
testStat, thresh = hsic_gam(X.reshape(len(X), 1), Y.reshape(len(Y), 1))
testStat / thresh
```

```
13.706487127525389
```

ノイズの大きさを少しずつ増やしていったときに、ピアソンの相関係数がどう変化するか計算してみましょう。

```Python
D = []
C = []
for d in range(100):
    X = np.linspace(-1, 1, 51)
    Y = X + np.random.rand(51,) * d / 10
    D.append(d)
    C.append(np.corrcoef(X, Y)[0, 1])

plt.plot(D, C)
```

![out14.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/0ac0a288-3205-5fdd-32b4-6ad0664bf8c0.png)

同様に、ノイズの大きさを少しずつ増やしていったときに、MICがどう変化するか計算してみましょう。

```Python

D = []
C = []
for d in range(100):
    X = np.linspace(-1, 1, 51)
    Y = X + np.random.rand(51,) * d / 10
    D.append(d)
    mine.compute_score(X, Y)
    C.append(mine.mic())

plt.plot(D, C)
```

![out15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/f20cdcc5-84ef-0978-0d76-27dbcdb1b0d0.png)

同様に、ノイズの大きさを少しずつ増やしていったときに、HSICがどう変化するか計算してみましょう。

```Python
D = []
C = []
for d in range(100):
    X = np.linspace(-1, 1, 51)
    Y = X + np.random.rand(51,) * d / 10
    D.append(d)
    testStat, thresh = hsic_gam(X.reshape(len(X), 1), Y.reshape(len(Y), 1))
    C.append(testStat / thresh)

plt.plot(D, C)
```

![out16.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/ad04c8c5-ea8f-1294-bfb7-381d218cf0fc.png)

# 二次関数
次は二次関数です。

```Python
import numpy as np

X = np.linspace(-1, 1, 51)
Y = X**2
plt.scatter(X, Y)
```

![out17.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/4d0ec064-d3df-2f31-533e-2e501535fe7d.png)

ピアソンの相関係数は

```Python
np.corrcoef(X, Y)[0, 1]
```

```
-2.3862043230836297e-17
```

MICは

```Python
mine.compute_score(X, Y)
mine.mic()
```

```
0.9997226475394071
```

HSICは

```Python
testStat, thresh = hsic_gam(X.reshape(len(X), 1), Y.reshape(len(Y), 1))
testStat / thresh
```

```
9.285886928766184
```

ノイズを加えてみましょう。

```Python
X = np.linspace(-1, 1, 51)
Y = X**2 + np.random.rand(51,)
plt.scatter(X, Y)
```

![out21.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/25a0d0e2-94a4-5bdc-0d86-40c8d3862e65.png)

ピアソンの相関係数は

```Python
np.corrcoef(X, Y)[0, 1]
```

```
0.04693076570744622
```

MICは

```Python
mine.compute_score(X, Y)
mine.mic()
```

```
0.5071787519579662
```

HSICは

```Python
testStat, thresh = hsic_gam(X.reshape(len(X), 1), Y.reshape(len(Y), 1))
testStat / thresh
```

```
5.247036645612692
```

ノイズの大きさを少しずつ増やしていったときに、ピアソンの相関係数がどう変化するか計算してみましょう。

```Python

D = []
C = []
for d in range(100):
    X = np.linspace(-1, 1, 51)
    Y = X**2 + np.random.rand(51,) * d / 10
    D.append(d)
    C.append(np.corrcoef(X, Y)[0, 1])

plt.plot(D, C)
```

![out25.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/ef061ecd-33a6-1d19-ebf2-b0cf46308171.png)

同様に、ノイズの大きさを少しずつ増やしていったときに、MICがどう変化するか計算してみましょう。

```Python

D = []
C = []
for d in range(100):
    X = np.linspace(-1, 1, 51)
    Y = X**2 + np.random.rand(51,) * d / 10
    D.append(d)
    mine.compute_score(X, Y)
    C.append(mine.mic())

plt.plot(D, C)
```

![out26.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/a1c4702e-8cc3-739b-6766-f0355434bd52.png)

同様に、ノイズの大きさを少しずつ増やしていったときに、HSICがどう変化するか計算してみましょう。

```Python
D = []
C = []
for d in range(100):
    X = np.linspace(-1, 1, 51)
    Y = X**2 + np.random.rand(51,) * d / 10
    D.append(d)
    testStat, thresh = hsic_gam(X.reshape(len(X), 1), Y.reshape(len(Y), 1))
    C.append(testStat / thresh)

plt.plot(D, C)
```

![out27.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/fc43f451-1001-846e-b8f6-cfe7c230a3a8.png)

ピアソンの相関係数が二次関数の関係性を検出できていないのに対して、MICもHSICも検出できていることがわかります。また、MICは小さなノイズにはほとんど影響を受けていないだろうことも分かります。

# 三角関数
三角関数も見てみましょう。

```Python
import numpy as np

X = np.linspace(-1, 1, 51)
Y = np.sin((X)*np.pi*2)
plt.scatter(X, Y)
```

![out39.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/94c14a90-4c5f-41a5-f46e-351a82b6929a.png)

ピアソンの相関係数は

```Python
np.corrcoef(X, Y)[0, 1]
```

```
-0.37651692742033543
```

MICは

```Python
mine.compute_score(X, Y)
mine.mic()
```

```
0.9997226475394071
```

HSICは

```Python
testStat, thresh = hsic_gam(X.reshape(len(X), 1), Y.reshape(len(Y), 1))
testStat / thresh
```

```
3.212604441715373
```

ノイズを加えてみましょう。

```Python
import numpy as np

X = np.linspace(-1, 1, 51)
Y = np.sin((X+0.25)*np.pi*2) + np.random.rand(51,)
plt.scatter(X, Y)
```

![out43.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/7fe9eb9f-c2b7-d3e4-c28a-fe497df50e37.png)

ピアソンの相関係数は

```Python
np.corrcoef(X, Y)[0, 1]
```

```
0.049378342546505014
```

MICは

```Python
mine.compute_score(X, Y)
mine.mic()
```

```
0.9997226475394071
```

HSICは

```Python
testStat, thresh = hsic_gam(X.reshape(len(X), 1), Y.reshape(len(Y), 1))
testStat / thresh
```

```
1.3135420547013614
```

ノイズをだんだん大きくしていったときの影響を調べます。ピアソンの相関係数は

```Python
D = []
C = []
for d in range(100):
    X = np.linspace(-1, 1, 51)
    Y = np.sin((X+0.25)*np.pi*2) + np.random.rand(51,) * d / 10
    D.append(d)
    C.append(np.corrcoef(X, Y)[0, 1])

plt.plot(D, C)
```

![out47.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/f8c989ca-d015-882a-dc34-e867351668de.png)

MICは

```Python
D = []
C = []
for d in range(100):
    X = np.linspace(-1, 1, 51)
    Y = np.sin((X+0.25)*np.pi*2) + np.random.rand(51,) * d / 10
    D.append(d)
    mine.compute_score(X, Y)
    C.append(mine.mic())

plt.plot(D, C)
```

![out48.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/8617b78c-3cec-460c-2b3f-f9a3dcf128da.png)

HSICは

```Python
D = []
C = []
for d in range(100):
    X = np.linspace(-1, 1, 51)
    Y = np.sin((X+0.25)*np.pi*2) + np.random.rand(51,) * d / 10
    D.append(d)
    testStat, thresh = hsic_gam(X.reshape(len(X), 1), Y.reshape(len(Y), 1))
    C.append(testStat / thresh)

plt.plot(D, C)
```

![out49.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/04e6243a-6555-54c2-b946-c4bc01e1f296.png)

ピアソンの相関係数が、ノイズ付きの三角関数の関係性を検出できなかったのは想定内でしたが、HSICも脱落してしまった感じですね。その中でもMICが関係性を検出できていたのは驚きでした。

# シグモイド関数
最後に、シグモイド関数です。

```Python
import numpy as np

X = np.linspace(-1, 1, 51)
Y = 1 / (1 + np.exp(-X*10))
plt.scatter(X, Y)
```

![out58.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/4e011540-f508-585f-544c-0249766a03d5.png)

ピアソンの相関係数は

```Python
np.corrcoef(X, Y)[0, 1]
```

```
0.9354020629807919
```

MICは

```Python
mine.compute_score(X, Y)
mine.mic()
```

```
0.9999999999999998
```

HSICは

```Python
testStat, thresh = hsic_gam(X.reshape(len(X), 1), Y.reshape(len(Y), 1))
testStat / thresh
```

```
27.433932271008874
```

ノイズを加えてみましょう。

```Python
import numpy as np

X = np.linspace(-1, 1, 51)
Y = 1 / (1 + np.exp(-X*10)) + np.random.rand(51,)
plt.scatter(X, Y)
```

![out62.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/b20f717a-1ce7-8764-4fe2-da052c739d0e.png)

ピアソンの相関係数は

```Python
np.corrcoef(X, Y)[0, 1]
```

```
0.7994507640548881
```

MICは

```Python
mine.compute_score(X, Y)
mine.mic()
```

```
0.9997226475394071
```

HSICは

```Python
testStat, thresh = hsic_gam(X.reshape(len(X), 1), Y.reshape(len(Y), 1))
testStat / thresh
```

```
12.218008194240534
```

ノイズを増やしていった影響を見てみましょう。ピアソンの相関係数は

```Python

D = []
C = []
for d in range(100):
    X = np.linspace(-1, 1, 51)
    Y = 1 / (1 + np.exp(-X*10)) + np.random.rand(51,) * d / 10
    D.append(d)
    C.append(np.corrcoef(X, Y)[0, 1])

plt.plot(D, C)
```

![out66.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/2257a5b1-349f-c4a3-61e7-7d2cdd7160c2.png)

MICは


```Python
D = []
C = []
for d in range(100):
    X = np.linspace(-1, 1, 51)
    Y = 1 / (1 + np.exp(-X*10)) + np.random.rand(51,) * d / 10
    D.append(d)
    mine.compute_score(X, Y)
    C.append(mine.mic())

plt.plot(D, C)
```

![out67.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/6b0862c2-3d08-ee12-fdfc-be7d04b2f823.png)

HSICは

```Python

D = []
C = []
for d in range(100):
    X = np.linspace(-1, 1, 51)
    Y = 1 / (1 + np.exp(-X*10)) + np.random.rand(51,) * d / 10
    D.append(d)
    testStat, thresh = hsic_gam(X.reshape(len(X), 1), Y.reshape(len(Y), 1))
    C.append(testStat / thresh)

plt.plot(D, C)
```

![out68.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/c0914caf-f09d-c1dc-4db5-f2cfb77a25e5.png)

いずれの指標も、関係性を検出できていますが、ノイズに対する強さで言えばMICが最強という感じです。

# 個人の感想です

MIC MIC にしてやんよ！