◇ 平均値の区間推定－ ロットアウトを防ごう－

標本標準偏差は、母標準偏差の 不偏推定量ではないが、母標準偏差の推定は、 近似的に標本標準偏差で行うことが多い ^{4 ) 5 )} 。

import numpy as np
import matplotlib.pyplot as plt
import random

sm = 52.2 # 平均（母平均）
ss = 9.5 # 標準偏差（母標準偏差）
sn = 10000 # 母数
en = 5 # 標本数
x = np.random.normal(loc=sm, scale=ss, size=sn)
sampled = random.sample(x.tolist(), en) #無作為抽出

fig = plt.figure()

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

ax1.hist(x)
ax2.hist(sampled)
plt.show()

average1 = np.mean(x)
stdev1 = np.std(x)

average2 = np.mean(sampled)
stdev2 = np.std(sampled)

print('inf',sm,ss)
print(sn,average1,stdev1)
print(en,average2,stdev2)

# ■■■ 平均値の有意差検定 ■■■

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 = 5 # 標本数
x = np.random.normal(loc=sm, scale=ss, size=sn)
sampled = random.sample(x.tolist(), en) #無作為抽出
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) #無作為抽出
ym = [ -0.5 for p in sampledm]

t, p = stats.ttest_ind(sampled, sampledm, alternative='two-sided')
#alternative=two-sided(両側検定), less(小なり片側検定), greater(大なり片側検定)
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")
ax1.axvline(np.mean(sampled),c="r")
ax1.text(np.mean(sampled), sn/5, " $\mu_\mathrm{f}$ = " + "{:.1f}".format(sm))

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, " $\overline{x}_\mathrm{f}$ = " + "{:.1f}".format(np.mean(sampled)))

ax2.hist(xm,color="cyan")
ax2.axvline(np.mean(sampledm),c="b")
ax2.text(np.mean(sampledm), sn/5, " $\mu_\mathrm{m}$ = " + "{:.1f}".format(smm))

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, " $\overline{x}_\mathrm{m}$ = " + "{:.1f}".format(np.mean(sampledm)))


ax3.text(30, 0.15, "$t$ = {:.3f}".format(t))
ax3.text(30, -0.15, "$p$ = {:.3f}".format(p))

plt.show()


# ■■■ 平均値の有意差検定 ■■■