variance.py

# coding: utf-8
# Your code here!
import numpy as np
data=np.array([[4.87,4.93,4.86,4.85],[4.86,4.90,4.85,4.86],[4.90,4.89,4.85,4.84],[4.87,4.91,4.81,4.86],[4.85,4.92,4.83,4.89]])
# data = np.array([[33,31,33],[30,29,31],[33,28,32],[29,29,32],[32,27,36]])
data_flatten = data.flatten()
data_mean = data_flatten.mean()
print("全体平均:{}".format(data_mean))

#列の平均
col_mean=[]
for i in range(len(data[0])):
    print(data[:,i])
    col_mean.append(data[:,i].mean())
print("列の平均:{}".format(col_mean))
  
 #列の効果
col_effect=[]
for i in range(len(col_mean)):
    col_effect.append(col_mean[i]-data_mean)
print("列の効果:{}".format(col_effect))

#誤差の計算
error_col=[]
error=[]
for i in range(len(data[0])):
    for j in range(len(data[:,i])):
        error_tmp = data[:,i][j] - col_mean[i]
        error_col.append(error_tmp)
    error.append(error_col)
    error_col=[]
print("各列ごとの誤差:{}".format(error))

#効果の分散を求める
sq_wa = 0
for i in range(len(col_effect)):
    sq_wa = sq_wa + len(data)*(col_effect[i]**2)
print(sq_wa)
dof = len(data[0])-1
col_effect_variance = sq_wa / dof
print("効果の分散:{}".format(col_effect_variance))

#誤差の分散を求める
error = np.array(error)
error = error.flatten()
sq_wa=0
for i in range(len(error)):
    sq_wa = sq_wa + error[i]**2
dof = len(error.flatten()) - len(data[0])
error_variance = sq_wa/dof
print("誤差の分散:{}".format(error_variance))

#分散比
variance_ration = col_effect_variance / error_variance
print("分散比:{}".format(variance_ration))

#F検定
#今回は自由度 (2,16)
#http://ktsc.cafe.coocan.jp/distributiontable.pdf