🐉

相関係数ハック

2023/03/14に公開

統計

numpy

tech

狙った相関係数になるように新しくデータを作るって方法あるんだっけ？

ってことで相関係数ハックしてみる。

考え方

相関係数といえば、これ

r = \frac{ \frac{1}{n} \Sigma_{i} (x_i - \bar{x})(y_i - \bar{y})} { \sqrt{\frac{1}{n} \Sigma_{i} (x_i - \bar{x})^2} \sqrt{\frac{1}{n} \Sigma_{i} (y_i - \bar{y})^2} }

今回は、これにガウスノイズを加えていくとする。
ガウスノイズを

n_i \sim N(0, \sigma)

とすると、

y_i = x_i + n_i

となる。　さらに、

\bar{y} = \bar{x}

でもある。

相関係数の式に代入すると、

r = \frac{ \frac{1}{n} \Sigma_{i} (x_i - \bar{x})(x_i + n_i - \bar{x})} { \sqrt{\frac{1}{n} \Sigma_{i} (x_i - \bar{x})^2} \sqrt{\frac{1}{n} \Sigma_{i} (x_i + n_i - \bar{x})^2} }

これをnumpy風にかくと、

r = \frac{ np.mean( (xs - xs.mean()) * (xs + ns - xs.mean())) }{ np.std(xs) * np.std(xs + ns) }

xs : input

ns : noise

となる。（これを書いててysのままでも良い気もしてきたが、そこはツッコミなしで）

なので、このrが狙った相関係数となるようなsの値を探索すれば良いのではと思った。

つまり、アルゴリズムとしては、↓のような感じになればいいんじゃないかと思った。

## 大域的なsの探索
for s in [search_min, search_max]:
    r_tmp = {ni ~ N(0, s) で上式から、rを計算}
    best_s = {r_tmpが探索範囲内で最小の時のs}
   
## sを固定した上で、何回かランダムにサンプリングして相関係数を探す
## ただし、狙った総関係数にならない場合もある（r_diff >= epsilon）　ので、最大の試行回数を予め決めておく
while (r_diff < epsilon) and (iter_n < max_iter):
    ns ~ N(0, best_s)
    ys = xs + ns
    r_tmp = corr(xs, ys)
    r_diff = (r_tmp - r_target)
    iter_n ++;

pythonコード

ということで、雑にpythonで書いてみる。（ちょっと変数名が違うが、そこは許してください。）

tqdmで処理の様子を可視化していたりするので、コードは上のpseudoコードとは少し違う。
最初のscaleの探索の上限下限はなんとなく、この範囲ならいけるっしょという勘（）

import numpy as np
import matplotlib.pyplot as plt
from tqdm.auto import tqdm

np.random.seed(42)

def get_corr_data(xs, corr=0.8, search_n=1000, eps=1e-6, iter_n=10_000, verbose=0):

    _data_n = len(xs)

    ## global scale search
    search_range = np.linspace(0.001 * np.ptp(xs), np.ptp(xs), search_n)
    best_diff = np.inf
    scale = -1

    pbar = tqdm(search_range, total=len(search_range), disable=verbose==0)

    for _s in pbar:
        ns = np.random.normal(0, _s, _data_n)
        _corr = np.mean( (xs - xs.mean()) * (xs + ns - xs.mean())) / np.std(xs) / np.std(xs + ns)
        _diff = abs(corr - _corr)

        if _diff < best_diff:
            best_diff = _diff
            scale = _s

        pbar.set_description(f'scale search : {best_diff = }  {scale = }')

    best_diff = np.inf
    ys_new = None

    ## correlation search
    pbar =  tqdm(range(iter_n), total=iter_n, disable=verbose==0)

    for iter in pbar:
        ns = np.random.normal(0, scale, _data_n)
        _ys = xs.copy() + ns
        _corr = np.corrcoef(xs, _ys)[0, 1]
        _diff = abs(corr - _corr)

        if _diff < best_diff:
            best_diff = _diff
            best_corr = _corr
            ys = _ys.copy()

        pbar.set_description(f'corr search : {best_corr = }')

        if best_diff < eps:
            break

    return ys

適当なxsで試す

データ数10,000で適当なsin波で試してみる。

data_n = 10_000
freq = 2

xs = np.sin(2 * np.pi * np.arange(data_n) / data_n * freq)

ちょっと可視化も含めて関数化しておく。

def sim1(xs):
    ys = get_corr_data(xs, corr=0.8, iter_n=10_000, eps=1e-5, search_n=1000, verbose=1)

    corr = np.corrcoef(xs, ys)[0, 1]
    print(f'{corr =}')

    fig, ax = plt.subplots(1, 2, figsize=(20, 5))
    ax[0].scatter(xs, ys)
    ax[0].set_title('xs vs ys')

    ax[1].plot(np.arange(data_n), xs, alpha=0.5, label='xs')
    ax[1].plot(np.arange(data_n), ys, alpha=0.5, label='ys')
    ax[1].legend()
    plt.show()

そして、これで試してみると、

sim1(xs)

うまくいってるんじゃないかな？？

他にもいろいろなxsを用意して試してみる。

for xs in [
    np.random.normal(0.4, 0.1, data_n),
    np.sin(2 * np.pi * np.arange(data_n) / data_n * freq) + np.linspace(-1, 1, data_n),
    np.sin(2 * np.pi * np.arange(data_n) / data_n * freq) * np.exp(-np.linspace(0, 2, data_n))
]:
    sim1(xs)

ガウスノイズもOKぽい。

にょきにょき上がっていくsin波もOKぽい。

高校の数学でよくみるこれもうまくいってるぽい。

まとめ

ということで、相関係数ハックできたと信じてる。

考え方

pythonコード

適当なxsで試す

まとめ

Discussion