Up

第５回　確率モデルと検定――前編――

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

リスト１　平均と標準偏差

import numpy as np

A = np.array([168, 171, 168, 167, 164, 172, 154, 174, 166, 166])

B = np.array([185, 181, 180, 187, 176, 182, 171, 176, 176, 184])

meanA = A.mean()

meanB = B.mean()

sdA = A.std()

sdB = B.std()

print('A国：　平均 = ', meanA, '　　SD = ', sdA)

print('B国：　平均 = ', meanB, ' SD = ', sdB)

リスト２　サンプリングデータの変動

import numpy as np

import matplotlib.pyplot as plt

import scipy.stats

import seaborn

mu = 170

sgm = 6

data = scipy.stats.norm.rvs(loc = mu, scale = sgm, size = 10)

print(data)

data = []

for i in range(5):

data.append(scipy.stats.norm.rvs(loc = mu, scale = sgm, size = 10))

Values = []

Groups = []

for i in range(5):

Values += [v for v in data[i]]

Groups += ['S-{}'.format(i + 1)] * len(data[i])

seaborn.stripplot(x = 'G', y = 'V', data = {'G': Groups, 'V': Values})

plt.show()

seaborn.boxplot(x = 'G', y = 'V', data = {'G': Groups, 'V': Values})

plt.show()

リスト３　t検定

import numpy as np

A = np.array([168, 171, 168, 167, 164, 172, 154, 174, 166, 166])

B = np.array([185, 181, 180, 187, 176, 182, 171, 176, 176, 184])

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

t, df = calc_t(A, B)

print('t = ',t, '　　自由度 = ', df)

リスト４　帰無仮説の下でのｔ分布のシミュレーション

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

mu = 170

sgm = 6

M = 10

N_Exp = 10000

t_values = []

for t in range(N_Exp):

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

Y = scipy.stats.norm.rvs(loc = mu, 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, bins = t_range, 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.show()

Up