比率の比較分析（ｐ値とベイズ分析）
２×２クロス表の分析
実行時エラーと対応（2024.09.03）
Ubuntu 22/WSLでの実行時エラー：
「cannot import name ’gaussian’ from ‘scipy.signal’
」というエラーが表示された。これは、scipyをバージョン1.12.0のものに置き換えると解決した。
置き換えは、インストールと同じコマンド
conda install
scipy==1.12.0
で行った。
Ubuntu 24/WSLでの実行時エラー：
Scipyに関するエラーは、Ubuntu
22の場合と同様、scipy 1.12.0をインストールすることで対応できた。
その後、実行すると
「Could not load the Qt platform plugin “xcb”」なるエラーが表示されるようになった。これは、表示デバイスがらみのエラーのようであるので、Ubuntu 24/WSLでは、StanによるMCMCサンプリングまでとして、後の処理はWindowsで行うことにした。例を、このウェブサイトの後半に示す。
Spiegelhalter(2019)は、腕組み（arm crossing）において左腕が上になる確率に性差があるかどうかの検定例を取り上げている。
女性14名において5名が左腕が上であり、男性40名において17名が左腕が上であった。カイ2乗検定とFisher’s Exact検定によるｐ値が扱われている。カイ２乗検定とFisher’s Exact検定のPythonスクリプトを書いてみた。また、これらのランダマイゼーション検定のスクリプトも試みた。これらは、ｐ値を求めるものであるが、ｐ値は、帰無仮説に対して算出されるデータの確率である。これに対して、ベイズ分析ではデータに対してモデル（パラメータ）の確率が求められる。ベイズ分析のスクリプトも試みた。カイ２乗検定の統計量は交互作用のみが対象となっているが、ベイズ分析では主効果も扱うことができる。
スクリプトファイルは圧縮ファイルc10files.zipにまとめた。
ライブラリScipyの関数scipy.stats.chi2_contingencyを用いるスクリプト例をリスト１のように作成した。
import numpy as np
import scipy.stats as ss
#
# Spiegelhalter(2019), Table 10.2
# Female Male
data = [[5, 17], # Left arm on top
[9, 23]] # Right arm on top
result = ss.chi2_contingency(data).pvalue
print('p-value =',result)
リスト１のスクリプトを実行すると、次の出力を得る。
p-value =
0.8975643879419304
Spiegelhalter(2019)のp.270に、P-valueが0.90であると記されている値が確認できた。カイ２乗検定の統計量は離散値であるが、連続量の分布関数であるカイ２乗分布による近似が用いられている。
次にFisherのexact検定を試みる。
Fisher’s exact検定は、Conover (1999),pp.188-189に従ってスクリプトを作成した。
|
表１ ２×２分割表 |
|
|
|
|
|
列１ |
列２ |
和 |
|
行１ |
|
|
|
|
行２ |
|
|
|
|
和 |
|
N-c |
|
表１の分割表において、左上のセルの値が
である確率が超幾何分布で与えられる。
ここで、の値が小さいとき、表１の行２列２の値が負になることがあるが、超幾何分布を与える関数は、このとき値0を返すので、の変域は、上式で考えてプログラミングができる。
Spiegelhalter(2019)のデータTable 10.2に対してFisher’s exact検定を行うプログラム例をリスト２に示す。
Fisher’s exact検定を行うスクリプトは、このウェブサイトにも示した。
import numpy as np
import scipy.stats as ss
#
# Spiegelhalter(2019), Table
10.2
#
Female Male
data = [[5, 17], # Left arm on top
[9,
23]] # Right arm on top
data = np.array(data)
#
# Conover(1999),
pp.188-189.
#
r = np.sum(data, axis=1)[0]
c = np.sum(data, axis=0)[0]
N = np.sum(data)
t_obs = data[0,0]
Ps = []
# List of probability
masses
for x in range(np.min([r,c])+1):
Ps.append(ss.hypergeom.pmf(x, N, r, c))
cumPs = np.cumsum(Ps)
Lower_tailed = cumPs[t_obs]
Upper_tailed = 1.0 - cumPs[t_obs-1] if
t_obs > 0 else 1.0
Two_tailed = 2 * np.min([Lower_tailed,
Upper_tailed])
print('Lower_tailed =', Lower_tailed)
print('Upper_tailed =', Upper_tailed)
print('Two_tailed =', Two_tailed)
リスト２のスクリプトを実行すると、以下の出力を得る。
Lower_tailed = 0.4528027112448938
Upper_tailed = 0.774794279497039
Two_tailed = 0.9056054224897876
両側検定におけるｐ値は、約0.91で、上のカイ2乗検定の値と大差ない。
カイ２乗分布による近似がどのようであるかを調べるために、ランダマイゼーションによるシミュレーションを行ってみた。
カイ2乗分布による近似の様子を調べるために、ランダマイゼーション検定を行ってみた。比較のために、Fisher’s exact検定もランダマイゼーション検定で行った。スクリプトをリスト３に示す。
リスト３のスクリプトを実行すると、Fisher’s exact検定の統計量、表１の、のヒストグラムが表示される。シミュレーションの試行数は、100000回である。
図R1
図R1のWindowを閉じると、図R2のグラフが表示される。
図R2
丸が理論値、超幾何分布の累積確率、を表し、菱形はランダマイゼーションにおけるFisher’s exact検定の統計量の累積密度を表している。極めてよく一致している。ランダマイゼーションはうまく行われたことを示していると考えられる。
図R2のウィンドウを閉じると、図R3のグラフが表示される。
図R3
統計量
のヒストグラムである。
図R3のウィンドウを閉じると、図R4の2つのグラフが表示される。
図R4
統計量の累積比率のグラフとカイ２乗分布の累積確率密度のグラフである。右側のグラフは、左側のグラフの左下部分を拡大表示したものである。大体近似できていると言えるが、グラフ上部に表示されているランダマイゼーションによるｐ値は0.759である。Fisher’s
exact検定の値0.9とはかなり離れている。
図Ｒ４のウィンドウを閉じると、プログラムは終了する。端末には以下の出力が表示されている。
p-value = 0.75874
p_chi2 = 0.65651214288276
p_chi2_Yates = 0.8975643879419304
p-valueは、ランダマイゼーション検定による値である。統計量の分布によるｐ値は0.9より小さい。
p_chi2は、カイ２乗分布の累積分布関数から求めたｐ値である。これも、0.9より小さい。
p_chi2_Yatesは、Yatesの修正を行った統計量のカイ2乗分布の累積分布関数から求めたｐである。実質的に0.9であるといえる。
Yatesの修正の分布をランダマイゼーションにより調べるスクリプトをリスト４のように用意した。
リスト４のスクリプトを実行すると、Yatesの修正のヒストグラムが表示される（図R5）。
図R5
図R5のウィンドウを閉じると、図R6のグラフが表示される。
図R6
修正の累積密度分布と累積カイ２乗分布のグラフである。グラフの上部に表示されているｐ値は、入力データSpiegelhalter(2019)のTable 10.2に対する修正とランダマイゼーション検定におけるｐ値である。ｐ値が1.000であるのは、Table 10.2に対する修正が最小値を取っていることを示すものである。
図R6のウィンドウを閉じると、プログラムは終了して、以下の出力が端末に表示される。
p-value(Yates' correction) by
simulation = 1.0
p_chi2_Yates = 0.8975643879419304
統計量修正に対応するカイ２乗分布によるｐ値は0.9で、Fisher’s exact検定によるｐ値に一致していると言える。しかし、修正の分布は、図Ｒ６に示されているようにカイ２乗による近似はよいとは言えない。また、ランダマイゼーション検定による修正のｐ値は１であることも注意すべき点である。
ｐ値は、帰無仮説に対して算出されるデータの確率であるが、これに対して、データに対してモデル（パラメータ）の確率を求める分析法としてベイズ分析がある。次に、ベイズ分析を試みる。
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
#
# Spiegelhalter(2019), Table
10.2
# Female Male
data = [[5, 17], # Left arm on top
[9,
23]] # Right arm on top
data = np.array(data)
n_f = np.sum(data, axis=0)[0]
n_m = np.sum(data, axis=0)[1]
n_L = np.sum(data, axis=1)[0]
n_R = np.sum(data, axis=1)[1]
N = np.sum(data)
#
# Expected values
#
e11 = n_f * n_L / N
e12 = n_m * n_L / N
e21 = n_f * n_R / N
e22 = n_m * n_R / N
#
# Statistic for
chi-squared test
#
chi2 = (data[0][0] - e11)**2 / e11 + \
(data[0][1] - e12)**2 / e12 + \
(data[1][0] - e21)**2 / e21 + \
(data[1][1] - e22)**2 / e22
def Yates(v):
return np.abs(v) - 0.5
chi2_Yates = (Yates(data[0][0] -
e11))**2 / e11 + \
(Yates(data[0][1] - e12))**2 / e12 + \
(Yates(data[1][0] - e21))**2 / e21 + \
(Yates(data[1][1] - e22))**2 / e22
n_simu = 100_000 # Number of simulation trials
results = []
# List of chi squared
values
resultsFisher = [] # List of statistics for Fisher's
exact test
# List of data for
randomization test
L_data = [1]*data[0][0] +
[0]*data[1][0] + [1]*data[0][1] + [0]*data[1][1]
L_data = np.array(L_data)
obs11_stat = np.sum(L_data[:n_f]) # Statistic for
Fisher's exact test
rng = np.random.default_rng()
#
# Simulation of randomization test
#
for i in range(n_simu):
rng.shuffle(L_data) # Shuffling for randomization test
obs11 = np.sum(L_data[:n_f])
obs12 =
np.sum(L_data[n_f-1:])
obs21 = n_f - obs11
obs22 = n_m - obs12
resultsFisher.append(obs11)
chi2_v = (obs11 - e11)**2 /
e11 + \
(obs12 - e12)**2 / e12 + \
(obs21 - e21)**2 / e21 + \
(obs22 - e22)**2 / e22
results.append(chi2_v)
plt.title("Distribution for
Fisher's Exact Test")
plt.hist(resultsFisher)
plt.show()
resultsFisher =
np.array(resultsFisher)
idxF = np.argsort(resultsFisher)
s_resultsFisher = resultsFisher[idxF]
subU_idxF = s_resultsFisher >=
obs11_stat # Upper tail
pvFisherU = np.mean(subU_idxF)
print('Fisher-Upper tailed =',
pvFisherU)
subL_idxF = s_resultsFisher <=
obs11_stat # Lower tail
pvFisherL = np.mean(subL_idxF)
print('Fisher-Lower tailed =',
pvFisherL)
pvFisherTwo = 2*pvFisherU if pvFisherU
< pvFisherL else 2*pvFisherL
print('Fisher-Two tailed =', pvFisherTwo)
r = np.sum(data, axis=1)[0]
c = np.sum(data, axis=0)[0]
N = np.sum(data)
FExact = []
for x in range(np.min([r, c])+1):
FExact.append(ss.hypergeom.cdf(x, N, r, c))
cum_F = np.cumsum(s_resultsFisher)
#
# Cumulative
distribution of samples
#
Fpsimu = []
for x in range(np.min([r,c])+1):
Fpsimu.append(np.mean(resultsFisher <= x))
plt.title("Fisher's exact
test")
plt.plot(range(np.min([r, c])+1),
Fpsimu, 'D', lw=0,
markersize=15, fillstyle='none', label='Simulation')
plt.plot(range(np.min([r, c])+1),
FExact, 's', lw=0,
markersize=15, fillstyle='none', label='Theory')
plt.legend()
plt.show()
plt.hist(results)
plt.title('Distribution for chi
squared test')
plt.xlabel('$\\chi^2$', fontsize=14)
plt.show()
plt.figure(figsize=(9,5))
plt.subplot(121)
ck = results >= (chi2 - 1.0e-7)
p_value = np.mean(ck)
# p-value
print('p-value =', p_value)
results = np.array(results)
idx = np.argsort(results) # Sorting for drawing cumulative
curve
s_results = results[idx]
# Drawing
cumulative distribution curve
plt.plot(s_results,
(np.arange(len(idx))+1)/len(idx),
label='Simulation')
# Cumulative
distribution of chi-squared distribution
chi2_c = ss.chi2.cdf(s_results, 1)
plt.plot(s_results, chi2_c,
label='$\\chi^2$')
plt.legend()
plt.title("Cumulative
distribution\n for chi-squared test")
plt.subplot(122)
sub_idx = s_results < 3*chi2
s_sub_results = s_results[sub_idx]
sub_chi2 = chi2_c[sub_idx]
plt.plot(s_sub_results, (np.arange(len(s_sub_results))+1)
/ len(idx),
label='Simulation')
plt.plot(s_sub_results, sub_chi2,
label='$\\chi^2$')
plt.title('Cumulative distribution\n
for chi squared test')
plt.legend(loc='lower right')
plt.suptitle(f"p-value =
{p_value:.3f}", fontsize=16)
plt.tight_layout()
plt.show()
p_chi2 = 1.0 - ss.chi2.cdf(chi2, 1)
print('p_chi2 =', p_chi2)
p_chi2_Yates = 1.0 -
ss.chi2.cdf(chi2_Yates, 1)
print('p_chi2_Yates =', p_chi2_Yates)
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
#
# Spiegelhalter(2019), Table
10.2
# Female Male
data = [[5, 17], # Left arm on top
[9,
23]] # Right arm on top
data = np.array(data)
n_f = np.sum(data, axis=0)[0]
n_m = np.sum(data, axis=0)[1]
n_L = np.sum(data, axis=1)[0]
n_R = np.sum(data, axis=1)[1]
N = np.sum(data)
#
# Expected values
#
e11 = n_f * n_L / N
e12 = n_m * n_L / N
e21 = n_f * n_R / N
e22 = n_m * n_R / N
#
#
Yates' Correction
#
def Yates(v):
return np.abs(v) - 0.5
chi2_Yates = (Yates(data[0][0] -
e11))**2 / e11 + \
(Yates(data[0][1] - e12))**2 / e12 + \
(Yates(data[1][0] - e21))**2 / e21 + \
(Yates(data[1][1] - e22))**2 / e22
n_simu = 100_000 # Number of simulation trials
results = []
# List of chi squared
values with Yates' correction
resultsFisher = [] # List of statistics for Fisher's
exact test
## List of data for
randomization test
L_data = [1]*data[0][0] +
[0]*data[1][0] + [1]*data[0][1] + [0]*data[1][1]
L_data = np.array(L_data)
rng = np.random.default_rng()
#
# Simulation of
randomization test
#
for i in range(n_simu):
rng.shuffle(L_data) # Shuffling for randomization test
obs11 = np.sum(L_data[:n_f])
obs12 =
np.sum(L_data[n_f-1:])
obs21 = n_f - obs11
obs22 = n_m - obs12
chi2_v = (Yates(obs11 -
e11))**2 / e11 + \
(Yates(obs12 - e12))**2 / e12 + \
(Yates(obs21 - e21))**2 / e21 + \
(Yates(obs22 - e22))**2 / e22
results.append(chi2_v)
plt.hist(results)
plt.title("Distribution for chi
squared test\nwith Yates' correction")
plt.xlabel('$\\chi^2$', fontsize=14)
plt.show()
plt.figure(figsize=(9,5))
plt.subplot(121)
ck = results >= (chi2_Yates -
1.0e-7)
p_value = np.mean(ck)
# p-value
print("p-value(Yates' correction)
by simulation =", p_value)
results = np.array(results)
idx = np.argsort(results) # Sorting for drawing cumulative
curve
s_results = results[idx]
# Drawing
cumulative distribution curve
plt.plot(s_results,
(np.arange(len(idx))+1)/len(idx),
label='Simulation')
# Cumulative
distribution of chi-squared distribution
chi2_c = ss.chi2.cdf(s_results, 1)
plt.plot(s_results, chi2_c,
label='$\\chi^2$')
plt.legend()
plt.title("Cumulative
distribution\n" +
"for chi-squared test with Yates' correction")
plt.subplot(122)
sub_idx = s_results < 0.6
s_sub_results = s_results[sub_idx]
sub_chi2 = chi2_c[sub_idx]
plt.plot(s_sub_results,
(np.arange(len(s_sub_results))+1) / len(idx),
label='Simulation')
plt.plot(s_sub_results, sub_chi2,
label='$\\chi^2$')
plt.title("Cumulative
distribution\n" +
"for chi squared test with Yates' corrction")
plt.legend(loc='lower right')
plt.suptitle(f"p-value =
{p_value:.3f}, $\\chi^2$(Yates) =
{chi2_Yates:.2f}",
fontsize=16)
plt.tight_layout()
plt.show()
p_chi2_Yates = 1.0 -
ss.chi2.cdf(chi2_Yates, 1)
print('p_chi2_Yates =', p_chi2_Yates)
2つの条件、列1と列２、のそれぞれにおいて2つの属性、行１と行２、のデータが表B1に示す個数得られたとする。
|
表B1 ２×２分割表 |
||
|
|
列１ |
列２ |
|
行１ |
|
|
|
行２ |
|
|
まず、通常のカイ２乗検定に合わせたモデルを設定する。その後、主効果の分析も含めたモデルによる分析を行う。
列１の条件において各データが属性行１である確率を
、列２の条件において各データが属性行１である確率を
として、次の２項分布を設定する。
上のモデルによるベイズ分析を行うStanスクリプトをリスト５のように用意した。
リスト５のStanスクリプトファイルxtable2by2.stanを利用してベイズ分析を行うPythonスクリプトをリスト６のように用意した。
リスト６のスクリプトを実行すると、MCMCサンプリングの終了後、図B1のトレース図が表示される。
図B1
図B1のWindowを閉じると、図B2のグラフが表示される。
図B2
列１と列２の確率との事後分布はかなり重なっており、差が無いようである。
図B2のWindowを閉じると、図B3のグラフが表示される。
図B3
確率の差
の事後分布のグラフである。差の左側の領域の確率が36.2％、右側の領域の確率が63.7％であるので、確率の差の事後分布のほぼ中央に値０が位置している。差が０でないとは言い難い事後分布である。
リスト５ ２×２分割表のベイズ分析のためのStanスクリプト（xtable2by2.stan）
data {
/*
o11 o12
o21 o22
*/
int o11;
int o12;
int o21;
int o22;
}
parameters {
real<lower=0, upper=1>
p1; // o11/(o11 + o21)
real<lower=0, upper=1>
p2; // o12/(o12 + o22)
}
model {
o11 ~ binomial(o11+o21, p1);
o12 ~ binomial(o12+o22, p2);
}
generated quantities {
real diff_p;
diff_p = p2 - p1;
}
import numpy as np
import stan
import pandas as pd
import arviz as az
import matplotlib.pyplot as plt
import seaborn as sb
#
# Spiegelhalter(2019), Table
10.2
#
Female Male
data = [[5, 17], # Left arm on top
[9,
23]] # Right arm on top
Data = {'o11':data[0][0], 'o12':data[0][1],
'o21':data[1][0], 'o22':data[1][1]}
with open('xtable2by2.stan') as
f_stan:
sm =
stan.build(f_stan.read(), data=Data)
fit = sm.sample()
infdata = az.from_pystan(posterior =
fit)
summary = az.summary(infdata)
print('\nSummary...\n', summary)
az.plot_trace(infdata)
plt.tight_layout()
plt.show()
dframe = fit.to_frame()
sb.kdeplot(dframe['p1'], color='b',
label='p1')
sb.kdeplot(dframe['p2'], color='g',
label='p2')
plt.xlabel('p1, p2', fontsize=14)
plt.legend(fontsize=16)
plt.tight_layout()
plt.show()
diff_p = np.array(dframe['diff_p'])
L_diff_p = diff_p <= 0.0
print(np.mean(L_diff_p))
R_diff_p = diff_p >= 0.0
print(np.mean(R_diff_p))
az.plot_kde(dframe['diff_p'])
plt.vlines([0.0], 0.0, 3.5)
plt.hlines([0.0], -0.4, 0.4)
plt.text(-0.152, 0.21,
f'{np.mean(L_diff_p)*100:.1f}%',
va='bottom', fontsize=16)
plt.text(0.05, 0.21,
f'{np.mean(R_diff_p)*100:.1f}%',
va='bottom', fontsize=16)
plt.xlabel('p2$-$p1', fontsize=16)
plt.tight_layout()
plt.show()
2つの条件、女性と男性、のそれぞれにおいて2つの属性、左腕が上と右腕が上、のデータが表B２に示す個数得られたとする。主効果を考えやすくするために、Spiegelhalter(2019), Table 10.2の記述に合わせた。
|
表B２ ２×２分割表（表B1の改変） |
||
|
|
女性 |
男性 |
|
左腕が上 |
|
|
|
右腕が上 |
|
|
まず、通常のカイ２乗検定に合わせたモデルは以下のように設定できる。
女性の条件において各データが左腕が上である確率を、男性において左腕が上である確率をとして、次の２項分布を設定する。
いま、表B2の回答者（受講生）が女性である確率を
で表し、２項分布を適用する。すなわち、
とおく。
上記のモデルによるベイズ分析を行うStanスクリプトをリスト７のように用意した。
リスト７のStanスクリプトを用いてベイズ分析を行うスクリプトをリスト８のように用意した。
リスト８のスクリプトを実行すると、MCMCサンプリング後に、まずトレース図が表示される（図M1）。
図M1
図M1のWindowを閉じると、図M2のグラフが表示される。
図M2
２つの確率の事後分布は、かなり重なっている。
図M2のWindowsを閉じると、図M3のグラフが表示される。
図M3
２つの確率の差、
、の分布はほぼ中央に０を含んでいる。
図M3のWindowを閉じると、図M4のグラフが表示される。
図M4
性差の主効果を表す確率の事後分布である。分布の領域が0.5以下の範囲に収まっている。統計学的に女性の確率は0.5より小さいと言える分布である。
標準のカイ2乗検定は交互作用のみが対象であるが、分割表に合わせた確率モデルを設定したベイズ分析では主効果も扱うことができる。
リスト７ 主効果を考量したStanスクリプト（xtable_fvsm.stan）
data {
/*
o11 o12
o21 o22
*/
int o11;
int o12;
int o21;
int o22;
}
parameters {
real<lower=0, upper=1>
pfm; // (o11+o21)/(o11+o21+o12+o22)
real<lower=0, upper=1>
p1; // o11/(o11 + o21)
real<lower=0, upper=1>
p2; // o12/(o12 + o22)
}
model {
o11 + o21 ~
binomial(o11+o21+o12+o22, pfm);
o11 ~ binomial(o11+o21, p1);
o12 ~ binomial(o12+o22, p2);
}
generated quantities {
real diff_p;
diff_p = p2 - p1;
}
リスト８ リスト７のStanスクリプト（xtable_fvsm.stan）を用いたベイズ分析スクリプト
import numpy as np
import stan
import pandas as pd
import arviz as az
import matplotlib.pyplot as plt
import seaborn as sb
#
# Spiegelhalter(2019), Table
10.2
# Female Male
data = [[5, 17], # Left arm on top
[9,
23]] # Right arm on top
Data = {'o11':data[0][0],
'o12':data[0][1],
'o21':data[1][0], 'o22':data[1][1]}
with open('xtable_fvsm.stan') as
f_stan:
sm = stan.build(f_stan.read(),
data=Data)
fit = sm.sample()
infdata = az.from_pystan(posterior =
fit)
summary = az.summary(infdata)
print('\nSummary...\n', summary)
az.plot_trace(infdata)
plt.tight_layout()
plt.show()
dframe = fit.to_frame()
sb.kdeplot(dframe['p1'], color='b',
label='p1')
sb.kdeplot(dframe['p2'], color='g',
label='p2')
plt.xlabel('p1, p2', fontsize=14)
plt.legend(fontsize=16)
plt.tight_layout()
plt.show()
diff_p = np.array(dframe['diff_p'])
L_diff_p = diff_p <= 0.0
print('L_diff_p =',np.mean(L_diff_p))
R_diff_p = diff_p >= 0.0
print('R_diff_p =', np.mean(R_diff_p))
az.plot_kde(dframe['diff_p'])
plt.vlines([0.0], 0.0, 3.5)
plt.text(0.01, 3.3, 'p2$-$p1$=0$',
fontsize=16)
plt.hlines([0.0], -0.4, 0.4)
plt.text(-0.152, 0.21, f'{np.mean(L_diff_p)*100:.1f}%',
va='bottom', fontsize=16)
plt.text(0.05, 0.21,
f'{np.mean(R_diff_p)*100:.1f}%',
va='bottom', fontsize=16)
plt.xlabel('p2$-$p1', fontsize=16)
plt.tight_layout()
plt.show()
pfm = np.array(dframe['pfm'])
L_pfm = pfm <= 0.5
print('L_pfm =', np.mean(L_pfm))
R_pfm = pfm >= 0.5
print('R_pfm =', np.mean(R_pfm))
az.plot_kde(dframe['pfm'])
plt.vlines([0.5], 0.0, 3.5)
plt.hlines([0.0], 0.0, 1.0)
plt.text(0.2, 0.21,
f'{np.mean(L_pfm)*100:.1f}%',
va='bottom', fontsize=16)
plt.text(0.51, 0.21,
f'{np.mean(R_pfm)*100:.1f}%',
va='bottom', fontsize=16)
plt.text(0.5, 3.5, 'pfm=0.5',
fontsize=16)
plt.xticks(np.arange(0, 1.01, 0.1))
plt.xlabel('pfm', fontsize=16)
plt.tight_layout()
plt.show()
Conover, W. J. (1999). Practical nonparametric statistics, third edition. John Wiley & Sons, Inc.
岡本安晴（2019） いまさら聞けないPythonでデータ分析．丸善出版
Spiegelhalter, D. (2019). The art of statistics: Learning from data. Pelican
Ubuntu 24/WSLでQtがらみのエラーが出る（2024.9.3）ので、Ubuntu 24/WSLでのStanによるベイズ分析は、MCMCサンプリングまでとして、結果のグラフ表示はWindowsで行った。例を以下に示す。
上の「ベイズ分析」のスクリプトをリストUW1とリストUW2に分割した例、および「主効果の分析を含めた場合」のスクリプトをリストUW3とリストUW4に分割した例を作成してみた。
リストUW1のスクリプトをUbuntu/WSLで実行してから、リストUW2のスクリプトをWindowsで実行すると、図UW!から始まって、グラフが順番に表示される。
図UW1
リストUW3のスクリプトをUbuntu/WSLで実行してから、リストUW4のスクリプトをWindowsで実行すると、図UW2のグラフから順番に表示される。
図UW2
import numpy as np
import stan
import pandas as pd
import arviz as az
#import matplotlib.pyplot as plt
#import seaborn as sb
import pickle
#
# Spiegelhalter(2019), Table
10.2
# Female Male
data = [[5, 17], # Left arm on top
[9,
23]] # Right arm on top
Data = {'o11':data[0][0],
'o12':data[0][1],
'o21':data[1][0], 'o22':data[1][1]}
with open('xtable2by2.stan') as
f_stan:
sm =
stan.build(f_stan.read(), data=Data)
fit = sm.sample()
infdata = az.from_pystan(posterior =
fit)
summary = az.summary(infdata)
print('\nSummary...\n', summary)
with open('infdata.pkl', 'wb') as f:
pickle.dump(infdata,f)
#az.plot_trace(infdata)
#plt.tight_layout()
#plt.show()
dframe = fit.to_frame()
with open('dframe.pkl', 'wb') as f:
pickle.dump(dframe, f)
import numpy as np
#import stan
import pandas as pd
import arviz as az
import matplotlib.pyplot as plt
import seaborn as sb
import pickle
#
# Spiegelhalter(2019), Table
10.2
#
Female Male
data = [[5, 17], # Left arm on top
[9,
23]] # Right arm on top
with open('infdata.pkl', 'rb') as f:
infdata = pickle.load(f)
az.plot_trace(infdata)
plt.tight_layout()
plt.show()
with open('dframe.pkl', 'rb') as f:
dframe = pickle.load(f)
sb.kdeplot(dframe['p1'], color='b',
label='p1')
sb.kdeplot(dframe['p2'], color='g',
label='p2')
plt.xlabel('p1, p2', fontsize=14)
plt.legend(fontsize=16)
plt.tight_layout()
plt.show()
diff_p = np.array(dframe['diff_p'])
L_diff_p = diff_p <= 0.0
print(np.mean(L_diff_p))
R_diff_p = diff_p >= 0.0
print(np.mean(R_diff_p))
az.plot_kde(dframe['diff_p'])
plt.vlines([0.0], 0.0, 3.5)
plt.hlines([0.0], -0.4, 0.4)
plt.text(-0.152, 0.21,
f'{np.mean(L_diff_p)*100:.1f}%',
va='bottom', fontsize=16)
plt.text(0.05, 0.21,
f'{np.mean(R_diff_p)*100:.1f}%',
va='bottom', fontsize=16)
plt.xlabel('p2$-$p1', fontsize=16)
plt.tight_layout()
plt.show()
import numpy as np
import stan
import pandas as pd
import arviz as az
#import matplotlib.pyplot as plt
#import seaborn as sb
import pickle
#
# Spiegelhalter(2019), Table
10.2
#
Female Male
data = [[5, 17], # Left arm on top
[9,
23]] # Right arm on top
Data = {'o11':data[0][0],
'o12':data[0][1],
'o21':data[1][0], 'o22':data[1][1]}
with open('xtable_fvsm.stan') as
f_stan:
sm =
stan.build(f_stan.read(), data=Data)
fit = sm.sample()
infdata = az.from_pystan(posterior =
fit)
summary = az.summary(infdata)
print('\nSummary...\n', summary)
with open('infdata_f.pkl', 'wb') as f:
pickle.dump(infdata, f)
#az.plot_trace(infdata)
#plt.tight_layout()
#plt.show()
dframe = fit.to_frame()
with open('dframe_f.pkl', 'wb') as f:
pickle.dump(dframe, f)
import numpy as np
#import stan
import pandas as pd
import arviz as az
import matplotlib.pyplot as plt
import seaborn as sb
import pickle
#
# Spiegelhalter(2019), Table
10.2
# Female Male
data = [[5, 17], # Left arm on top
[9,
23]] # Right arm on top
with open('infdata_f.pkl', 'rb') as f:
infdata = pickle.load(f)
az.plot_trace(infdata)
plt.tight_layout()
plt.show()
with open('dframe_f.pkl', 'rb') as f:
dframe = pickle.load(f)
sb.kdeplot(dframe['p1'], color='b',
label='p1')
sb.kdeplot(dframe['p2'], color='g',
label='p2')
plt.xlabel('p1, p2', fontsize=14)
plt.legend(fontsize=16)
plt.tight_layout()
plt.show()
diff_p = np.array(dframe['diff_p'])
L_diff_p = diff_p <= 0.0
print('L_diff_p =',np.mean(L_diff_p))
R_diff_p = diff_p >= 0.0
print('R_diff_p =', np.mean(R_diff_p))
az.plot_kde(dframe['diff_p'])
plt.vlines([0.0], 0.0, 3.5)
plt.text(0.01, 3.3, 'p2$-$p1$=0$', fontsize=16)
plt.hlines([0.0], -0.4, 0.4)
plt.text(-0.152, 0.21,
f'{np.mean(L_diff_p)*100:.1f}%',
va='bottom', fontsize=16)
plt.text(0.05, 0.21,
f'{np.mean(R_diff_p)*100:.1f}%',
va='bottom', fontsize=16)
plt.xlabel('p2$-$p1', fontsize=16)
plt.tight_layout()
plt.show()
pfm = np.array(dframe['pfm'])
L_pfm = pfm <= 0.5
print('L_pfm =', np.mean(L_pfm))
R_pfm = pfm >= 0.5
print('R_pfm =', np.mean(R_pfm))
az.plot_kde(dframe['pfm'])
plt.vlines([0.5], 0.0, 3.5)
plt.hlines([0.0], 0.0, 1.0)
plt.text(0.2, 0.21,
f'{np.mean(L_pfm)*100:.1f}%',
va='bottom', fontsize=16)
plt.text(0.51, 0.21,
f'{np.mean(R_pfm)*100:.1f}%',
va='bottom', fontsize=16)
plt.text(0.5, 3.5, 'pfm=0.5',
fontsize=16)
plt.xticks(np.arange(0, 1.01, 0.1))
plt.xlabel('pfm', fontsize=16)
plt.tight_layout()
plt.show()