分散の有意差検定をしよう

# ■■■ 母分散の有意差検定 ■■■

#帰無仮説H0:分散が等しい
#対立仮説H1:分散に有意差がある

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import random
from scipy import stats

sm = 52.2 # 平均（母平均）
ss = 9.5 # 標準偏差（母標準偏差）
sn = 10000 # 母数
en = 20 # 標本数
x = np.random.normal(loc=sm, scale=ss, size=sn)
sampled = random.sample(x.tolist(), en) #無作為抽出
vf = np.var(sampled,ddof=1) #ddof=1を指定して、標本分散ではなく不偏分散を求める
y = [ 0.5 for p in sampled]

smm = 64.8 # 平均（母平均）
ssm = 13.9 # 標準偏差（母標準偏差）
xm = np.random.normal(loc=smm, scale=ssm, size=sn)
sampledm = random.sample(xm.tolist(), en) #無作為抽出
vm = np.var(sampledm,ddof=1) #ddof=1を指定して、標本分散ではなく不偏分散を求める

ym = [ -0.5 for p in sampledm]

f = vm/vf if vm>vf else vf/vm #F値が1以上になるようF値を求める
dfn = en - 1 #自由度は標本数-1
dfd = en - 1 #自由度は標本数-1
p = 1 - stats.f.cdf(f, dfn=dfn, dfd=dfd) #累積分布関数の右側の面積
result = "Significant" if p<0.05 else "Not significant"

fig = plt.figure()

ax1 = fig.add_subplot(3, 1, 1)
ax2 = fig.add_subplot(3, 1, 2)
ax3 = fig.add_subplot(3, 1, 3)

ax1.hist(x, color="pink",alpha=0.2)
ax1.axvline(np.mean(x),c="r")
ax1.axvspan(xmin=np.mean(x)-np.std(x), xmax=np.mean(x)+np.std(x), color="pink", alpha=0.8)
ax1.text(sm, sn/5, " $\sigma_\mathrm{f}$ = " + "{:.1f}".format(ss))

e2 = patches.Ellipse(xy=(np.mean(sampled), +0.5), width=np.std(sampled)*6, height=0.2, fc='pink')
ax3.add_patch(e2)
ax3.axvspan(xmin=np.mean(sampled)-np.std(sampled), xmax=np.mean(sampled)+np.std(sampled), color="pink", alpha=0.3)
ax3.axvline(np.mean(sampled),c="r")
ax3.scatter(sampled,y,c="r")
ax3.text(np.mean(sampled), +0.25, " $s_\mathrm{f}$ = " + "{:.1f}".format(np.std(sampled,ddof=1)))

ax2.hist(xm,color="cyan",alpha=0.2)
ax2.axvline(np.mean(xm),c="b")
ax2.axvspan(xmin=np.mean(xm)-np.std(xm), xmax=np.mean(xm)+np.std(xm), color="cyan", alpha=0.8)
ax2.text(smm, sn/5, " $\sigma_\mathrm{m}$ = " + "{:.1f}".format(ssm))

e1 = patches.Ellipse(xy=(np.mean(sampledm), -0.5), width=np.std(sampledm)*6, height=0.2, fc='cyan')
ax3.add_patch(e1)
ax3.axvspan(xmin=np.mean(sampledm)-np.std(sampledm), xmax=np.mean(sampledm)+np.std(sampledm), color="cyan", alpha=0.3)
ax3.axvline(np.mean(sampledm),c="b")
ax3.scatter(sampledm,ym,c="b")
ax3.text(np.mean(sampledm), -0.35, " $s_\mathrm{m}$ = " + "{:.1f}".format(np.std(sampledm,ddof=1)))

ax3.text(30, 0, result)
ax3.text(30, 0.2, "$F$ = {:.3f}".format(f))
ax3.text(30, -0.2, "$p$ = {:.3f}".format(p))

plt.show()


# ■■■ 母分散の有意差検定 ■■■