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

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

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

MIC (Maximal Information Coefficient) を計算するPythonパッケージ

次のようにして pip install できます。

!pip install minepy

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

from minepy import MINE
mine = MINE()

HSIC (Hilbert-Schmidt Independence Criterion) を計算するPythonパッケージ

次のように git clone します。

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

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

from HSIC.HSIC import hsic_gam

一次関数

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

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

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

ピアソンの相関係数は

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

1.0

MICは

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

0.9999999999999998

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

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

22.610487654066624

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

import numpy as np

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

ピアソンの相関係数は

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

0.9018483832673769

MICは

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

0.9997226475394071

HSICは

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

13.706487127525389

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

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)

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

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)

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

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)

二次関数

次は二次関数です。

import numpy as np

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

ピアソンの相関係数は

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

-2.3862043230836297e-17

MICは

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

0.9997226475394071

HSICは

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

9.285886928766184

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

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

ピアソンの相関係数は

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

0.04693076570744622

MICは

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

0.5071787519579662

HSICは

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

5.247036645612692

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

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)

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

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)

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

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)

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

三角関数

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

import numpy as np

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

ピアソンの相関係数は

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

-0.37651692742033543

MICは

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

0.9997226475394071

HSICは

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

3.212604441715373

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

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)

ピアソンの相関係数は

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

0.049378342546505014

MICは

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

0.9997226475394071

HSICは

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

1.3135420547013614

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

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)

MICは

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)

HSICは

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)

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

シグモイド関数

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

import numpy as np

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

ピアソンの相関係数は

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

0.9354020629807919

MICは

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

0.9999999999999998

HSICは

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

27.433932271008874

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

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)

ピアソンの相関係数は

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

0.7994507640548881

MICは

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

0.9997226475394071

HSICは

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

12.218008194240534

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

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)

MICは

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)

HSICは

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)

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

個人の感想です

MIC MIC にしてやんよ！