t_test.py

#!python3
#-*- coding: utf-8 -*-
# v0.5 2014-07-23
# v0.6 2016-07-01


import numpy as np
#import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import math


# 有意水準５%値
yuui_5 =((0, 0),
					(1, 12.706),
					(2, 4.303),
					(3, 3.182),
					(4, 2.776),
					(5, 2.571),
					(6, 2.447),
					(7, 2.365),
					(8, 2.306),
					(9, 2.262),
					(10, 2.226),
					(11, 2.201),
					(12, 2.179),
					(13, 2.160),
					(14, 2.145),
					(15, 2.131),
					(16, 2.120),
					(17, 2.110),
					(18, 2.101),
					(19, 2.093),
					(20, 2.086),
					(21, 2.086),
					(22, 2.074),
					(23, 2.069),
					(24, 2.064),
					(25, 2.060),
					(26, 2.056),
					(27, 2.052),
					(28, 2.048),
					(29, 2.045),
					(30, 2.042),
					(40, 2.021),
					(60, 2.000),
					(120, 1.980),
					(1000, 1.960),
					)
					
					
# 簡易データ分析
def summary(data1, data2):
	print('----- summary -----')
	print('標本1= ', data1)
	print('標本数= ', len(data1))
	print('最小値= ', np.min(data1))
	print('最大値= ', np.max(data1))
	print('中央値= ', np.median(data1))
	print('平均= ', np.mean(data1))
	print('分散= ', np.var(data1))
	print('-----')
	print('標本2= ', data2)
	print('標本数= ', len(data2))
	print('最小値= ', np.min(data2))
	print('最大値= ', np.max(data2))
	print('中央値= ', np.median(data2))
	print('平均= ', np.mean(data2))
	print('分散= ', np.var(data2))


# 不偏分散
def fuhen_var(np_arr):
	mean = np.mean(np_arr)
	sum = 0
	for i in np_arr:
		sum = sum + math.pow(i-mean, 2)
		return sum / (len(np_arr)-1)

# 標本標準誤差
def hyoujyun_gosa(np_arr):
	mean = np.mean(np_arr)
	sum = 0
	for i in np_arr:
		sum = sum + math.pow(i-mean, 2)
		fuhen = sum / (len(np_arr)-1)
	return math.sqrt(fuhen / len(np_arr))

# 散布図表示
def plot(data1, data2):
	plt.plot(data1, 'bo')
	plt.plot(data2, 'gx')
	plt.margins(0.2)
	plt.show()
	plt.close()

# ヒストグラム表示
def hist(data1, data2):
	plt.hist(data1)
	plt.hist(data2)
	plt.show()
	plt.close()

# t検定
def t_test(data1, data2):
	bunsan = ((np.var(data1)*len(data1)) + (np.var(data2)*len(data2))) / ((len(data1)-1) + (len(data2)-1))
	gosa = math.sqrt(bunsan * ((1.0/len(data1)) + (1.0/len(data2))))
	heikin_sa = np.mean(data1)-np.mean(data2)
	#
	jiyuudo = len(data1)-1 + len(data2)-1
	if jiyuudo > 30 and jiyuudo < 40:
		jiyuudo = 30
	elif jiyuudo > 40 and jiyuudo < 60:
		jiyuudo = 40
	elif jiyuudo > 60 and jiyuudo < 120:
		jiyuudo = 60
	elif jiyuudo > 120:
		jiyuudo = 1000
	#	
	t = heikin_sa / gosa
	#
	t_5 = yuui_5[jiyuudo][1]
	#
	print('------- t検定 --------')
	print('自由度= ', jiyuudo)
	print('t = ', t)
	print('(t分布表)有意水準5%値 = ',t_5)
	print(' ')
	print('95%信頼区間 = ')
	print((heikin_sa - (t_5 * gosa)), '~', (heikin_sa + (t_5 * gosa)))


def main():
	data1 = np.array(DATA1)
	data2 = np.array(DATA2)
	summary(data1, data2)
	# 
	while True:
		print(' ')
		print('s)ummary, p)lot, h)istgram, t)_test, q)uit = ?')
		inp = input()
		if inp == 'q':
			print('****************')
			print('** 終了しました **')
			print('****************')
			break
		elif inp == 's':
			summary(data1, data2)
		elif inp == 'p':
			plot(data1, data2)
		elif inp == 'h':
			hist(data1, data2)
		elif inp == 't':
			t_test(data1, data2)



if __name__ == '__main__':
	# 対象データ
	#DATA1 = [70,75,70,85,90,70,80,75]
	#DATA2 = [85,80,95,70,80,75,80,90]
	#DATA1 = [80, 75, 80, 95, 90, 80, 85, 90]
	#DATA2 = [75,65,80,85,75,80,80,70]
	#DATA1 = [7,8,10,5,8,7,9,5,6,9,10,6,7,8,7,9,10,6]
	#DATA2 = [9,9,6,10,9,8,10,7,9,10,6,8,9,9,10,7,8,8,10,9]
	DATA1 = [15,18,17,16,16,18,12,17,14,16,16,13,16,15,15,15,12,14,16]
	DATA2 = [15,12,16,13,15,14,15,16,15,13,13,15,16,14,15,17,14,15,15]
	# 
	main()