# -*- coding: UTF-8 -*-
from math import *
import numpy as np
############################################
# Γ（ｘ）の計算（ガンマ関数，近似式）
# ier : =0 : normal
# =-1 : x=-n (n=0,1,2,･･･)
# return : 結果
# coded by Y.Suganuma
############################################
def gamma(x, ier) :
ier[0] = 0
if x > 5.0 :
v = 1.0 / x
s = ((((((-0.000592166437354 * v + 0.0000697281375837) * v + 0.00078403922172) * v - 0.000229472093621) * v - 0.00268132716049) * v + 0.00347222222222) * v + 0.0833333333333) * v + 1.0
g = 2.506628274631001 * exp(-x) * pow(x,x-0.5) * s
else :
err = 1.0e-20
w = x
t = 1.0
if x < 1.5 :
if x < err :
k = int(x)
y = float(k) - x
if abs(y) < err or abs(1.0-y) < err :
ier[0] = -1
if ier[0] == 0 :
while w < 1.5 :
t /= w
w += 1.0
else :
if w > 2.5 :
while w > 2.5 :
w -= 1.0
t *= w
w -= 2.0
g = (((((((0.0021385778 * w - 0.0034961289) * w + 0.0122995771) * w - 0.00012513767) * w + 0.0740648982) * w + 0.0815652323) * w + 0.411849671) * w + 0.422784604) * w + 0.999999926
g *= t
return g
################################################
# 標準正規分布N(0,1)の計算（P(X = x), P(X < x)）
# w : P(X = x)
# return : P(X < x)
################################################
def normal(x, w) :
# 確率密度関数（定義式）
w[0] = exp(-0.5 * x * x) / sqrt(2.0*pi)
# 確率分布関数（近似式を使用）
y = 0.70710678118654 * abs(x)
z = 1.0 + y * (0.0705230784 + y * (0.0422820123 + y * (0.0092705272 + y * (0.0001520143 + y * (0.0002765672 + y * 0.0000430638)))))
P = 1.0 - z ** (-16.0)
if x < 0.0 :
P = 0.5 - 0.5 * P
else :
P = 0.5 + 0.5 * P
return P
#######################################
# ｆ分布の計算（P(X = ff), P(X < ff)）
# dd : P(X = ff)
# df1,df2 : 自由度
# return : P(X < ff)
#######################################
def f(ff, dd, df1, df2) :
if ff < 1.0e-10 :
ff = 1.0e-10
x = ff * df1 / (ff * df1 + df2)
if df1%2 == 0 :
if df2%2 == 0 :
u = x * (1.0 - x)
pp = x
ia = 2
ib = 2
else :
u = 0.5 * x * sqrt(1.0-x)
pp = 1.0 - sqrt(1.0-x)
ia = 2
ib = 1
else :
if df2%2 == 0 :
u = 0.5 * sqrt(x) * (1.0 - x)
pp = sqrt(x)
ia = 1
ib = 2
else :
u = sqrt(x*(1.0-x)) / pi
pp = 1.0 - 2.0 * atan2(sqrt(1.0-x), sqrt(x)) / pi
ia = 1
ib = 1
if ia != df1 :
for i1 in range(ia, df1-1, 2) :
pp -= 2.0 * u / i1
u *= x * (i1 + ib) / i1
if ib != df2 :
for i1 in range(ib, df2-1, 2) :
pp += 2.0 * u / i1
u *= (1.0 - x) * (i1 + df1) / i1
dd[0] = u / ff
return pp
############################################
# 二分法による非線形方程式(f(x)=0)の解
# fn : f(x)を計算する関数名
# x1,x2 : 初期値
# eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）
# eps2 : 終了条件２（｜f(x(k))｜＜eps2）
# max : 最大試行回数
# ind : 実際の試行回数
# （負の時は解を得ることができなかった）
# return : 解
# coded by Y.Suganuma
############################################
def bisection(fn, x1, x2, eps1, eps2, max, ind) :
x0 = 0.0
f1 = fn(x1)
f2 = fn(x2)
if f1*f2 > 0.0 :
ind[0] = -1
else :
ind[0] = 0
if f1*f2 == 0.0 :
if f1 == 0.0 :
x0 = x1
else :
x0 = x2
else :
sw = 0
while sw == 0 and ind[0] >= 0 :
sw = 1
ind[0] += 1
x0 = 0.5 * (x1 + x2)
f0 = fn(x0)
if abs(f0) > eps2 :
if ind[0] <= max :
if abs(x1-x2) > eps1 and abs(x1-x2) > eps1*abs(x2) :
sw = 0
if f0*f1 < 0.0 :
x2 = x0
f2 = f0
else :
x1 = x0
f1 = f0
else :
ind[0] = -1
return x0
############################################
# Newton法による非線形方程式(f(x)=0)の解
# fn : f(x)を計算する関数名
# dfn : f(x)の微分を計算する関数名
# x0 : 初期値
# eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）
# eps2 : 終了条件２（｜f(x(k))｜＜eps2）
# max : 最大試行回数
# ind : 実際の試行回数
# （負の時は解を得ることができなかった）
# return : 解
# coded by Y.Suganuma
############################################
def newton(fn, dfn, x0, eps1, eps2, max, ind) :
x1 = x0
x = x1
ind[0] = 0
sw = 0
while sw == 0 and ind[0] >= 0 :
sw = 1
ind[0] += 1
g = fn(x1)
if abs(g) > eps2 :
if ind[0] <= max :
dg = dfn(x1)
if abs(dg) > eps2 :
x = x1 - g / dg
if abs(x-x1) > eps1 and abs(x-x1) > eps1*abs(x) :
x1 = x
sw = 0
else :
ind[0] = -1
else :
ind[0] = -1
return x
############################
# ｆ分布の計算
# coded by Y.Suganuma
############################
##########################################
# 1.0 - p - P(X>x)（関数値, 標準正規分布）
##########################################
def normal_f(x) :
y = np.empty(1, np.float)
return 1.0 - p - normal(x, y)
################################################################
# 標準正規分布N(0,1)のｐ％値（P(X > u) = 0.01p）（二分法を使用）
# ind : >= 0 : normal（収束回数）
# = -1 : 収束しなかった
################################################################
def p_normal(ind) :
u = bisection(normal_f, -7.0, 7.0, 1.0e-6, 1.0e-10, 100, ind)
return u
#######################################
# 1.0 - p - P(X > x)（関数値, f 分布）
#######################################
def f_f(x) :
y = np.empty(1, np.float)
return f(x, y, dof1, dof2) - 1.0 + p
#################################
# P(X = x)（関数の微分, f 分布）
#################################
def f_df(x) :
y = np.empty(1, np.float)
z = f(x, y, dof1, dof2)
return y[0]
######################################
# ｆ分布のｐ％値（P(X > u) = 0.01p）
# ind : >= 0 : normal（収束回数）
# = -1 : 収束しなかった
######################################
def p_f(ind) :
MAX = 340
ff = 0.0
sw = 0
# 初期値計算の準備
# while文は，大きな自由度によるガンマ関数の
# オーバーフローを避けるため
while sw >= 0 :
df1 = 0.5 * (dof1 - sw)
df2 = 0.5 * dof2
a = 2.0 / (9.0 * (dof1 - sw))
a1 = 1.0 - a
b = 2.0 / (9.0 * dof2)
b1 = 1.0 - b
yq = p_normal(ind)
e = b1 * b1 - b * yq * yq
if e > 0.8 or dof1+dof2-sw <= MAX :
sw = -1
else :
sw += 1
if (dof1-sw) == 0 :
sw = -2
if sw == -2 :
ind[0] = -1
else :
# f0 : 初期値
if e > 0.8 :
x = (a1 * b1 + yq * sqrt(a1*a1*b+a*e)) / e
f0 = x ** 3.0
else :
y1 = float(dof2) ** (df2-1.0)
y2 = float(dof1) ** df2
x = gamma(df1+df2, ind) / gamma(df1, ind) / gamma(df2, ind) * 2.0 * y1 / y2 / p
f0 = x ** (2.0/dof2)
# ニュートン法
ff = newton(f_f, f_df, f0, 1.0e-6, 1.0e-10, 100, ind)
return ff
#######################################
# 分散分析（一元配置法と二元配置法）
# method : 1 or 2
# Np[0] : 因子1の水準数
# [1] : 因子2の水準数
# N : データ数
# name[0] : 因子1の名前
# [1] : 因子2の名前
# a : a ％値
# x : データ
# coded by Y.Suganuma
#######################################
def aov(method, Np, N, name, a, x) :
global p, dof1, dof2
# 一元配置法
if method == 1 :
xi = np.empty(Np[0], np.float)
for i1 in range(0, Np[0]) :
xi[i1] = 0.0
for i2 in range(0, N) :
xi[i1] += x[i1][0][i2]
xi[i1] /= N
xa = 0.0
for i1 in range(0, Np[0]) :
for i2 in range(0, N) :
xa += x[i1][0][i2]
xa /= (Np[0] * N)
SP = 0.0
for i1 in range(0, Np[0]) :
SP += (xi[i1] - xa) * (xi[i1] - xa)
SP *= N
SE = 0.0
for i1 in range(0, Np[0]) :
for i2 in range(0, N) :
SE += (x[i1][0][i2] - xi[i1]) * (x[i1][0][i2] - xi[i1])
VP = SP / (Np[0] - 1)
VE = SE / (Np[0] * (N - 1))
FP = VP / VE
p = 0.01 * a
dof2 = Np[0] * (N - 1)
sw = np.empty(1, np.int)
print("全変動： 平方和 " + str(SP+SE) + " 自由度 " + str(Np[0]*N-1))
dof1 = Np[0] - 1
print(name[0] + "の水準間： 平方和 " + str(SP) + " 自由度 " + str(Np[0]-1) + " 不偏分散 " + str(VP) + " F " + str(FP) + " " + str(a) + "％値 " + str(p_f(sw)))
print("水準内： 平方和 " + str(SE) + " 自由度 " + str(Np[0]*(N-1)) + " 不偏分散 " + str(VE))
# 二元配置法
else :
xi = np.empty(Np[0], np.float)
for i1 in range(0, Np[0]) :
xi[i1] = 0.0
for i2 in range(0, Np[1]) :
for i3 in range(0, N) :
xi[i1] += x[i1][i2][i3]
xi[i1] /= (Np[1] * N)
xj = np.empty(Np[1], np.float)
for i1 in range(0, Np[1]) :
xj[i1] = 0.0
for i2 in range(0, Np[0]) :
for i3 in range(0, N) :
xj[i1] += x[i2][i1][i3]
xj[i1] /= (Np[0] * N)
xij = np.empty((Np[0], Np[1]), np.float)
for i1 in range(0, Np[0]) :
for i2 in range(0, Np[1]) :
xij[i1][i2] = 0.0
for i3 in range(0, N) :
xij[i1][i2] += x[i1][i2][i3]
xij[i1][i2] /= N
xa = 0.0
for i1 in range(0, Np[0]) :
for i2 in range(0, Np[1]) :
for i3 in range(0, N) :
xa += x[i1][i2][i3]
xa /= (Np[0] * Np[1] * N)
SP = 0.0
for i1 in range(0, Np[0]) :
SP += (xi[i1] - xa) * (xi[i1] - xa)
SP *= (Np[1] * N)
SQ = 0.0
for i1 in range(0, Np[1]) :
SQ += (xj[i1] - xa) * (xj[i1] - xa)
SQ *= (Np[0] * N)
SI = 0.0
for i1 in range(0, Np[0]) :
for i2 in range(0, Np[1]) :
SI += (xij[i1][i2] - xi[i1] - xj[i2] + xa) * (xij[i1][i2] - xi[i1] - xj[i2] + xa)
SI *= N
SE = 0.0
for i1 in range(0, Np[0]) :
for i2 in range(0, Np[1]) :
for i3 in range(0, N) :
SE += (x[i1][i2][i3] - xij[i1][i2]) * (x[i1][i2][i3] - xij[i1][i2])
VP = SP / (Np[0] - 1)
VQ = SQ / (Np[1] - 1)
VI = SI / ((Np[0] - 1) * (Np[1] - 1))
VE = SE / (Np[0] * Np[1] * (N - 1))
FP = VP / VE
FQ = VQ / VE
FI = VI / VE
p = 0.01 * a
dof2 = Np[0] * Np[1] * (N - 1)
sw = np.empty(1, np.int)
print("全変動： 平方和 " + str(SP+SQ+SI+SE) + " 自由度 " + str(Np[0]*Np[1]*N-1))
dof1 = Np[0] - 1
print(name[0] + "の水準間： 平方和 " + str(SP) + " 自由度 " + str(Np[0]-1) + " 不偏分散 " + str(VP) + " F " + str(FP) + " " + str(a) + "％値 " + str(p_f(sw)))
dof1 = Np[1] - 1
print(name[1] + "の水準間： 平方和 " + str(SQ) + " 自由度 " + str(Np[1]-1) + " 不偏分散 " + str(VQ) + " F " + str(FQ) + " " + str(a) + "％値 " + str(p_f(sw)))
dof1 = (Np[0] - 1) * (Np[1] - 1)
print("相互作用： 平方和 " + str(SI) + " 自由度 " + str((Np[0]-1)*(Np[1]-1)) + " 不偏分散 " + str(VI) + " F " + str(FI) + " " + str(a) + "％値 " + str(p_f(sw)))
print("水準内： 平方和 " + str(SE) + " 自由度 " + str(Np[0]*Np[1]*(N-1)) + " 不偏分散 " + str(VE))
p = 0.0
dof1 = 1
dof2 = 2
----------------------------------
# -*- coding: UTF-8 -*-
import numpy as np
import sys
from math import *
from function import aov, f, gamma, normal, bisection, newton
############################
# 分散分析
# coded by Y.Suganuma
############################
name = ["", ""]
Np = np.array([1, 1], np.int)
line = sys.stdin.readline()
ss = line.split()
method = int(ss[0]) # 因子の数
N = int(ss[1]) # 各水準におけるデータ数
a = float(ss[2]) # 有意水準（％）
if method == 1 or method == 2 :
for i1 in range(0, method) :
line = sys.stdin.readline()
ss = line.split()
name[i1] = ss[0]
Np[i1] = int(ss[1])
x = np.empty((Np[0], Np[1], N), np.float)
for i1 in range(0, Np[1]) :
for i2 in range(0, N) :
line = sys.stdin.readline()
ss = line.split()
for i3 in range(0, Np[0]) :
x[i3][i1][i2] = float(ss[i3])
aov(method, Np, N, name, a, x)
else :
print("一元配置法，または，二元配置法だけです！")
-------- 一元配置法に対するデータ例（コメント部分を除いて下さい）--------
1 6 5 // 因子の数 各水準におけるデータ数 有意水準（％）
工場 3 // 因子の名前 水準の数
3.1 4.7 5.1 // 各水準に対する1番目のデータ
4.1 5.6 3.7 // 各水準に対する2番目のデータ
3.3 4.3 4.5 // 各水準に対する3番目のデータ
3.9 5.9 6.0 // 各水準に対する4番目のデータ
3.7 6.1 3.9 // 各水準に対する5番目のデータ
2.4 4.2 5.4 // 各水準に対する6番目のデータ
-------- 二元配置法に対するデータ例（コメント部分を除いて下さい）--------
2 3 5 // 因子の数 各水準におけるデータ数 有意水準（％）
薬剤 5 // 1番目の因子の名前 その水準の数
品種 2 // 2番目の因子の名前 その水準の数
3 4 12 -4 -4 // 1番目の因子の各水準に対して，2番目の因子の水準1に対する1番目のデータ
8 -8 31 12 19 // 1番目の因子の各水準に対して，2番目の因子の水準1に対する2番目のデータ
7 -5 8 0 23 // 1番目の因子の各水準に対して，2番目の因子の水準1に対する3番目のデータ
8 -10 9 10 15 // 1番目の因子の各水準に対して，2番目の因子の水準2に対する1番目のデータ
-5 11 26 -1 13 // 1番目の因子の各水準に対して，2番目の因子の水準2に対する2番目のデータ
10 -6 13 -7 -6 // 1番目の因子の各水準に対して，2番目の因子の水準2に対する3番目のデータ