相関係数のサンプリング・シミュレーション
真の相関係数が０の場合
相関が無い（真の相関係数ρ＝０）の場合の、サンプルサイズが小さいときのサンプル相関係数ｒの分布について調べる。サンプルサイズNが５の場合を扱うが、これは質的研究などにおいてサンプルをランダムにとって相関係数ｒを求めた場合、データ数が少ないと相関があるかのような印象を受けることがあるので、注意が必要であることを示すものである。シミュレーションに用いたスクリプトファイルは、圧縮ファイルfiles.zipにまとめた。
真の相関係数ρが０であるとき、サンプルサイズN=5の相関係数ｒを次式でｔに変換した値
は、自由度
のｔ分布に従う（岡本、2022）。
式（１）をｒについて解くと式（２）を得る。
式（２）のグラフをリスト１のスクリプトで描くと、図１のようになる。
図１
「相関係数の真値
」の条件の下で
である。すなわち、相関が無い場合でも、データ数がN=5と少ないときは、相関係数のサンプル値はかなり広い範囲、−0.878〜0.878、を取り得る。質的研究において、数人のデータから変量間の関係の有無を論ずるときは注意が必要である。
リスト１ 相関係数ｒの分布（r_smpl_distri.py）
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
N = 5
ts = np.linspace(-6, 6, 1000)
t2s = ts**2
r2s = t2s / ((N-2) + t2s)
rs = r2s**0.5
tp0975 = ss.t.ppf(0.975, N-2)
def t_r(t, N):
return ((t**2) / ((N-2) +
(t**2))) ** 0.5
r0975 = t_r(tp0975, N)
plt.plot([-tp0975, tp0975], [r0975,
r0975], label = f'r = {r0975:.3f}')
plt.plot([tp0975, tp0975], [-r0975,
r0975], label = f't = {tp0975:.3}')
plt.plot([-tp0975, tp0975], [-r0975,
-r0975], label = f'r = {-r0975:.3f}')
plt.plot([-tp0975, -tp0975], [-r0975,
r0975], label = f't = {-tp0975:.3}')
plt.plot(ts, rs)
plt.plot(ts, -rs)
plt.plot([-6, 6], [0, 0], c = 'k')
plt.xlabel('$t$', fontsize = 16)
plt.ylabel('$r$', fontsize = 16)
plt.legend(loc = 'lower right')
plt.title('$N = ${0}'.format(N),
fontsize = 20)
plt.savefig('Fig_r_smpls.png')
plt.show()
「真の相関係数ρ＝０」のもとで、サンプルサイズN=5のデータに対して50回シミュレーションを繰り返したときのρの事後分布を求めると図２のようになった。
図２
真の相関係数ρが０なので、事後分布はｒ＝０の周りに集まっているが、ｒ＝±１の近くに分布している事後分布も認められる。
各事後分布において、正の領域の確率Posと負の領域の確率Negを求め、ヒストグラムを描くと図３のようになった。
図３
正の領域、あるいは負の領域の確率が０あるいは１に近い度数が無視できない頻度であることが解る。すなわち、「真の相関係数ρ＝０」の場合であっても、サンプルサイズがN=5と少ないときは、ベイズ分析の場合も事後分布による相関係数の推定値が０であると言えない頻度が無視できないことが分かる。
リスト２のスクリプトによりサンプルに対してベイズ分析をStanで行っている。Stanによる事後分布からのサンプリングに基づいて、リスト３のスクリプトで事後分布の分析を行っている。
StanはPython用のインターフェースPyStan 3を用いている。PyStan 3の使用については、このウェブサイトで説明している。
PyStan３をWindowsで用いるときは、WSL (Windows Subsystem for Linux)を利用する。WSLについては、このウェブサイトで説明している。
すなわち、Windowsにおいては、リスト２のスクリプトはAnaconda/Ubuntu/WSL/WindowsのPython 3.10仮想環境で実行して、その後、リスト３のスクリプトをWindows上で実行するという手順になる。
なお、PyStanについては拙著＜岡本安晴（2019）＞で説明している。
リスト２ 事後分布をStanにより求めるスクリプト（PostDistCorr.py）
import numpy as np
import scipy.stats as ss
#import matplotlib.pyplot as plt
import seaborn as sb
import stan
import pickle
import arviz as az
model_script = """
data {
int
N;
vector[2] XY[N];
}
parameters {
real
mu1;
real
mu2;
real<lower = 0> sgm1;
real<lower = 0> sgm2;
real<lower=-1, upper=1> rho;
}
transformed parameters {
vector[2] mus;
matrix[2,2] cov;
mus[1] = mu1;
mus[2] = mu2;
cov[1][1] = sgm1^2;
cov[2][2] = sgm2^2;
cov[1][2] = sgm1*sgm2*rho;
cov[2][1] = cov[1][2];
}
model {
mu1
~ normal(0,10);
mu2
~ normal(0,10);
sgm1
~ exponential(0.01);
sgm2
~ exponential(0.01);
rho
~ uniform(-1, 1);
//for (i in 1:N){
// XY[i] ~
multi_normal(mus, cov);
XY ~
multi_normal(mus, cov);
//}
}
"""
N = 5
rhos = []
for i in range(50):
smpls =
ss.multivariate_normal.rvs(mean= [0, 0],
cov = [[1, 0.0], [0.0, 1]], size = N)
Data = {'N':N, 'XY':smpls}
sm =
stan.build(model_script, data = Data)
fit = sm.sample()
infer_data =
az.from_pystan(posterior=fit, posterior_model=sm)
smry =
az.summary(infer_data)
print(smry)
rhos.append(fit['rho'])
with open('fits.pkl', 'wb') as f:
pickle.dump(rhos, f)
リスト３ リスト３の実行結果の描画スクリプト（DrawFit.py）
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import seaborn as sb
import pickle
import arviz as az
with open('fits.pkl', 'rb') as f:
fits = pickle.load(f)
for fit in fits:
sb.kdeplot(fit[0])
plt.xlabel('$r$', fontsize = 18)
plt.title('$N = ${0}'.format(5),
fontsize = 20)
plt.savefig('FigPosts.png')
plt.show()
p_pos = []
p_neg = []
for fit in fits:
cnt_pos = [1 if v > 0
else 0 for v in fit[0]]
cnt_neg = [1 if v < 0
else 0 for v in fit[0]]
p_pos.append(sum(cnt_pos) /
len(cnt_pos))
p_neg.append(sum(cnt_neg) /
len(cnt_neg))
plt.hist(p_pos, color = 'y', alpha =
0.5, label = 'Pos')
plt.hist(p_neg, color = 'b', alpha =
0.5, label = 'Neg')
plt.xlabel('Probability', fontsize =
16)
plt.legend()
plt.savefig('FigHists.png')
plt.show()
引用文献
岡本安晴（2022）データを使いこなすための統計入門．22世紀アート （電子書籍）
岡本安晴（2019）いまさら聞けないPythonでデータ分析．丸善出版