Up

第６回　確率モデルと検定――後編――

プログラム（Pythonスクリプト）の完全なリストを以下に上げる。スクリプトのファイルは圧縮ファイルfiles6.zipとしてまとめたので、クリックしてダウンロード解凍すればスクリプトファイルを開くことができる。プログラムの説明は連載本文を参照されたい。

リスト１　対立仮説の下でのｔ値のシミュレーション

import numpy as np

import matplotlib.pyplot as plt

import scipy.stats

def calc_t(X, Y):

M = len(X)

N = len(Y)

df = M + N - 2

meanX = X.mean()

meanY = Y.mean()

ss = (((X - meanX) ** 2).sum() + ((Y - meanY) ** 2).sum()) / df

ss_XY = ss * (M + N) / (M * N)

t = (meanX - meanY) / (ss_XY ** 0.5)

return t, df

muA = 170

muB = 180

sgm = 6

M = 10

N_Exp = 10000

t_values = []

for t in range(N_Exp):

X = scipy.stats.norm.rvs(loc = muA, scale = sgm, size = M)

Y = scipy.stats.norm.rvs(loc = muB, scale = sgm, size = M)

t, df = calc_t(X, Y)

t_values.append(t)

t_range = [v for v in np.arange(-4, 4, 0.25)]

p_values = [scipy.stats.t.pdf(t, 2 * M - 2) for t in t_range]

plt.hist(t_values, density = True)

plt.plot(t_range, p_values, linewidth = 3)

plt.plot([-5.45, -5.45], [0.0, 0.1], 'r-o', linewidth = 5, label = 't = -5.45')

plt.legend()

plt.title('$\mu$A = {0} $\mu$B = {1}'.format(muA, muB), fontsize = 16)

plt.show()

リスト２　棄却域の描画

import numpy as np

import matplotlib.pyplot as plt

import scipy.stats

df = 18

t_a2L = scipy.stats.t.ppf(0.025, df)

t_a2R = scipy.stats.t.ppf(0.975, df)

print('P(t<{0:.3f}) = 0.025, P(t>{1:.3f}) = 0.025'.format(t_a2L, t_a2R))

t_range = np.arange(-4.0, 4.0, 0.01)

p_values = scipy.stats.t.pdf(t_range, df)

plt.plot(t_range, p_values)

plt.hlines(0.0, -4.0, 4.0, color = 'k')

plt.vlines(t_a2L, 0.0, scipy.stats.t.pdf(t_a2L, df), color = 'r', linewidth = 2,

label = 'P(t<{0:<0.2f}) = 0.025'.format(t_a2L))

t_rangeL = np.arange(-4.0, t_a2L, 0.01)

plt.fill_between(t_rangeL, scipy.stats.t.pdf(t_rangeL, df), color = 'r', alpha = 0.5)

plt.vlines(t_a2R, 0.0, scipy.stats.t.pdf(t_a2R, df), color = 'm', linewidth = 2,

label = 'P(t>{0:<.2f}) = 0.025'.format(t_a2R))

t_rangeR = np.arange(t_a2R, 4.0, 0.01)

plt.fill_between(t_rangeR, scipy.stats.t.pdf(t_rangeR, df), color = 'm', alpha = 0.5)

plt.title('t distribution with df = 18', fontsize = 16)

plt.legend()

plt.show()

リスト３　ｔ値と確率

import scipy.stats

df = 18

t = -2.0

p = scipy.stats.t.cdf(t, df)

print('p = ', p)

Up