相関係数のサンプリング・シミュレーション
真の相関係数が0.3の場合
相関がある（真の相関係数ρ＝0.3）場合の、サンプルサイズが小さいときのサンプル相関係数ｒの分布について調べる。サンプルサイズNが５の場合を扱うが、これは質的研究などにおいてサンプルをランダムにとって相関係数ｒを求めた場合、データ数が少ないと相関が認められないことがあるので、注意が必要であることを示すものである。シミュレーションに用いたスクリプトファイルは、圧縮ファイルfiles0p3.zipにまとめた。
真の相関係数がρ＝0.3であるとき、サンプルサイズN=5のデータの相関係数を10000回のシミュレーションで求めると、図１のような分布になった。使用したスクリプトはリスト１に示すものである。
図１
ヒストグラムの横軸上に、下記の式（１）および（２）に基づいて求めた帰無仮説「ρ＝０」の採択域を示している。
真の相関係数ρが０であるとき、サンプルサイズN=5の相関係数ｒを次式でｔに変換した値
は、自由度
のｔ分布に従う（岡本、2022）。
式（１）をｒについて解くと式（２）を得る。
式（１）からｔ分布の棄却の臨界値を求めることができ、式（２）より対応するサンプル相関係数ｒの値が求まる。
サンプル相関係数の採択域は
と相関係数の値の範囲
とほぼ重なっていて、サンプル相関係数ｒの分布も採択域の上にほとんどがある。すなわち、真の相関係数がρ＝0.3であっても、帰無仮説「ρ＝０」が棄却されることはほとんどない。
リスト１ 相関係数ｒの分布（r_smpl_distri.py）
import numpy
as np
import scipy.stats
as ss
import matplotlib.pyplot
as plt
N = 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)
r_sample = []
for i in
range(10000):
smpls
= ss.multivariate_normal.rvs(mean= [0, 0],
cov = [[1, 0.3], [0.3, 1]], size = N)
results = ss.pearsonr(smpls.T[0], smpls.T[1])
r_sample.append(results[0])
plt.plot([-r0975, r0975], [0, 0], lw = 10, label = 'NonSig.Interval'
+
'\n[{0:.2f}, {1:.2f}]'.format(-r0975, r0975))
plt.hist(r_sample,
bins = 20)
plt.legend(fontsize =
16)
plt.xlabel('$r$', fontsize
= 20)
plt.title(r'$\rho=0.3,$' + ' ' + '$N=5$', fontsize
= 20)
plt.savefig('FigSmpls.png')
plt.show()
「真の相関係数ρ＝0.3」のもとで、サンプルサイズN=5のデータに対して50回のシミュレーションをリスト２のスクリプトにより繰り返し、得られた事後分布のサンプルをリスト３のスクリプトにより分析すると図２のようになった。
図２
真の相関係数ρ＝0.3に対応して、事後分布はｒの正の領域への偏りが認められるが、負の領域の分布も無視できない。
各事後分布において、正の領域の確率Posと負の領域の確率Negを求め、ヒストグラムを描くと図３のようになった。
図３
事後分布の負の領域の確率が0.05以下である頻度は少ないので、ρの値が正であると推定できることはあまり期待できない。
リスト２のスクリプトによりサンプルに対してベイズ分析を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);
XY ~
multi_normal(mus, cov);
}
"""
N = 5
rhos = []
for i in
range(50):
smpls
= ss.multivariate_normal.rvs(mean= [0, 0],
cov = [[1, 0.3], [0.3, 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でデータ分析．丸善出版