代码收藏家 方差分析代码(Python)
方差分析代码_Python
单因素方差分析
函数
import numpy as np
from scipy.stats import f,t
def One_way_ANOVA(dt,alpha,est):
"""
单因素方差分析,并输出单因素试验方差分析表
:param dt: (列表推导式) 行代表因素A的不同水平，每行代表该水平下的多次试验数据(元素个数可以不同)
:param alpha: 置信度
:param est: 需要参数估计的均值索引
:return: none
"""
r = len(dt) #因素A水平数
ni = np.zeros(r,dtype=int)
for i in range(r):
ni[i] = len(dt[i]) #因素A的不同水平Ai下的试验次数
n = np.sum(ni,dtype=int) #总试验次数
'将dt数据转化为numpy数组,数据缺失部分用0补上'
data = np.zeros((r,max(ni)))
for i,dt_row in enumerate(dt):
data[i, :len(dt_row)] = dt_row
'QA~误差平方和 QE~误差偏差平方和 QT~总偏差平方和'
Q = np.sum((np.sum(data,axis=1)**2) /ni)
P = (np.sum(data))**2 /n
R = np.sum(data**2)
QA = Q-P
QE = R-Q
QT = R-P
'fa~QA对应的自由度 fe~QE对应的自由度 ft~QT对应的自由度'
fa = r-1
fe = n-r
ft = n-1
'平均离差 + F检验量的值'
QA_bar = QA/fa
QE_bar = QE/fe
F = QA_bar/QE_bar
'画单因素试验方差分析表'
print(' 方差分析表')
print('=========================================================')
print("方差来源 平方和 自由度 均方 F比")
print('---------------------------------------------------------')
print('组间 ' + f'{QA:.2f}'+' '*(15-len(f'{QA:.2f}')) + f'{fa:d}'+' '*(12-len(f'{fa:d}'))
+ f'{QA_bar:.2f}'+' '*(12-len(f'{QA_bar:.2f}')) + f'{F:.2f}')
print('组内 ' + f'{QE:.2f}'+' '*(15-len(f'{QE:.2f}')) + f'{fe:d}'+' '*(12-len(f'{fe:d}'))
+ f'{QE_bar:.2f}'+' '*(12-len(f'{QE_bar:.2f}')))
print('总和 ' + f'{QT:.2f}'+' '*(15-len(f'{QT:.2f}')) + f'{ft:d}'+' '*(12-len(f'{ft:d}')))
print('=========================================================\n')
print("结论:")
'判断F是不是在拒绝域{F≥F_alpha(fa,fe)}里面'
F_alpha = f.isf(alpha,fa,fe) #F(fa,fe)的上α分位数
if F>F_alpha:
print(f"F={F:.2f} ≥ F_{alpha:.2f}({fa:d}, {fe:d})={F_alpha:.2f}，在显著性水平{alpha:.2f}上 拒绝 原假设H0。\n")
else:
print(f"F={F:.2f} < F_{alpha:.2f}({fa:d}, {fe:d})={F_alpha:.2f}，在显著性水平{alpha:.2f}上 接受 原假设H0。\n")
'参数估计'
if est is not None:
print("参数估计:")
Ai_m = np.sum(data,axis=1)/ni #单因素每个水平下的试验数据均值
r1 = est.shape[0] #需要进行参数估计的次数
for i in range(r1):
low = (Ai_m[est[i,0]-1]-Ai_m[est[i,1]-1]) - t.isf(alpha/2,n-r)*np.sqrt(QE_bar*(1/ni[est[i,0]-1]+1/ni[est[i,1]-1]))
upp = (Ai_m[est[i,0]-1]-Ai_m[est[i,1]-1]) + t.isf(alpha/2,n-r)*np.sqrt(QE_bar*(1/ni[est[i,0]-1]+1/ni[est[i,1]-1]))
print(f'μ{est[i,0]:d}-μ{est[i,1]:d}的置信水平为{1-alpha:.2f}的置信区间为: ({low:.2f},{upp:.2f})')
例4
if __name__=="__main__": #主函数
data1 = [[231.0,232.8,227.6,228.3,224.7,225.5,229.3,230.3],
[222.8,224.5,218.5,220.2],
[224.3,226.1,221.4,223.6]]
alpha = 0.01
est = None
One_way_ANOVA(data1, alpha, est)
例5
if __name__=="__main__": #主函数
data1 = [[256,222,280,298],
[244,300,290,275],
[250,277,230,322],
[288,280,315,259],
[206,212,220,212]]
alpha = 0.05
est = np.array([[1,5]])
One_way_ANOVA(data1, alpha, est)
两因素非重复方差分析
函数
import numpy as np
from scipy.stats import f,t
def Two_way_ANOVA(data,alpha):
"""
两因素非重复方差分析,输出试验方差分析表
:param data: np.array([r,s]) 行代表因素A的不同水平，列代表因素B的不同水平，位置ij的元素代表(Ai,Bj)下的单次试验数据
:param alpha: 置信度
:return: none
"""
r,s = data.shape #因素A的水平数，因素B的水平数
n = r*s #总试验次数
'QA QB QE QT'
P = (np.sum(data)**2) /n
R = np.sum(data**2)
Q1 = np.sum(np.sum(data,axis=1)**2) /s
Q2 = np.sum(np.sum(data,axis=0)**2) /r
QA = Q1-P
QB = Q2-P
QE = R-Q1-Q2+P
QT = R-P
'fa~QA对应的自由度 fb~QB对应的自由度 fe~QE对应的自由度 ft~QT对应的自由度'
fa = r-1
fb = s-1
fe = (r-1)*(s-1)
ft = r*s-1
'平均离差 + F检验量'
QA_bar = QA/fa
QB_bar = QB/fb
QE_bar = QE/fe
FA = QA_bar/QE_bar
FB = QB_bar/QE_bar
'画两因素非重复试验方差分析表'
print(' 方差分析表')
print('=========================================================')
print("方差来源 平方和 自由度 均方 F比")
print('---------------------------------------------------------')
print('因素A ' + f'{QA:.2f}'+' '*(15-len(f'{QA:.2f}')) + f'{fa:d}'+' '*(12-len(f'{fa:d}'))
+ f'{QA_bar:.2f}'+' '*(12-len(f'{QA_bar:.2f}')) + f'{FA:.2f}')
print('因素B ' + f'{QB:.2f}'+' '*(15-len(f'{QB:.2f}')) + f'{fb:d}'+' '*(12-len(f'{fb:d}'))
+ f'{QB_bar:.2f}'+' '*(12-len(f'{QB_bar:.2f}')) + f'{FB:.2f}')
print('误差 ' + f'{QE:.2f}'+' '*(15-len(f'{QE:.2f}')) + f'{fe:d}'+' '*(12-len(f'{fe:d}'))
+ f'{QE_bar:.2f}'+' '*(12-len(f'{QE_bar:.2f}')))
print('总和 ' + f'{QT:.2f}'+' '*(15-len(f'{QT:.2f}')) + f'{ft:d}'+' '*(12-len(f'{ft:d}')))
print('=========================================================\n')
print("结论:")
'判断FA,FB是不是在拒绝域{FA≥F_alpha(fa,fe)},{FB≥F_alpha(fb,fe)}里面'
FA_alpha = f.isf(alpha,fa,fe) # FA(fa,fe)的上α分位数
if FA>FA_alpha:
print(f"FA={FA:.2f} ≥ FA_{alpha:.2f}({fa:d}, {fe:d})={FA_alpha:.2f}，在显著性水平{alpha:.2f}上 拒绝 原假设H01。")
else:
print(f"FA={FA:.2f} < FA_{alpha:.2f}({fa:d}, {fe:d})={FA_alpha:.2f}，在显著性水平{alpha:.2f}上 接受 原假设H01。")
FB_alpha = f.isf(alpha,fb,fe) # FB(fb,fe)的上α分位数
if FB > FB_alpha:
print(f"FB={FB:.2f} ≥ FB_{alpha:.2f}({fb:d}, {fe:d})={FB_alpha:.2f}，在显著性水平{alpha:.2f}上 拒绝 原假设H02。")
else:
print(f"FB={FB:.2f} < FB_{alpha:.2f}({fb:d}, {fe:d})={FB_alpha:.2f}，在显著性水平{alpha:.2f}上 接受 原假设H02。")
例1
if __name__=="__main__": #主函数
data2 = np.array([[63.1,63.9,65.6,66.8],
[65.1,66.4,67.8,69.0],
[67.2,71.0,71.9,73.5]])
alpha = 0.01
Two_way_ANOVA(data2, alpha)
两因素等重复方差分析
函数
import numpy as np
from scipy.stats import f,t
def Two_way_ANOVA_repeated(data,alpha):
"""
两因素等重复方差实验分析,输出试验方差分析表
:param data: np.array([r,s,t]) 行代表因素A的不同水平，列代表因素B的不同水平，位置ijk的元素代表(Ai,Bj)下的第k次试验数据
:param alpha: 置信度
:return: none
"""
r,s,t = data.shape # 因素A的水平数，因素B的水平数，水平(Ai,Bj)下的试验次数
n = r*s*t # 总试验次数
data = np.transpose(data, (2, 0, 1)) #将np.array([r,s,t])转换成np.array([t,r,s])以方便后续计算
'QA QB QAB QE QT'
P = (np.sum(data)**2) /n
W = np.sum(data**2)
U = np.sum(np.sum(np.sum(data,axis=0),axis=1)**2) /(s*t)
V = np.sum(np.sum(np.sum(data,axis=0),axis=0)**2) /(r*t)
R = np.sum(np.sum(data,axis=0)**2) /t
QA = U-P
QB = V-P
QAB = R-U-V+P
QE = W-R
QT = W-P
'fa~QA对应的自由度 fb~QB对应的自由度 fab~QAB对应的自由度 fe~QE对应的自由度 ft~QT对应的自由度'
fa = r-1
fb = s-1
fab = (r-1)*(s-1)
fe = r*s*(t-1)
ft = r*s*t-1
'平均离差 + F检验量'
QA_bar = QA/fa
QB_bar = QB/fb
QAB_bar = QAB/fab
QE_bar = QE/fe
FA = QA_bar/QE_bar
FB = QB_bar/QE_bar
FAB = QAB_bar/QE_bar
'画两因素非重复试验方差分析表'
print(' 方差分析表')
print('=========================================================')
print("方差来源 平方和 自由度 均方 F比")
print('---------------------------------------------------------')
print('因素A ' + f'{QA:.2f}'+' '*(15-len(f'{QA:.2f}')) + f'{fa:d}'+' '*(12-len(f'{fa:d}'))
+ f'{QA_bar:.2f}'+' '*(12-len(f'{QA_bar:.2f}')) + f'{FA:.2f}')
print('因素B ' + f'{QB:.2f}'+' '*(15-len(f'{QB:.2f}')) + f'{fb:d}'+' '*(12-len(f'{fb:d}'))
+ f'{QB_bar:.2f}'+' '*(12-len(f'{QB_bar:.2f}')) + f'{FB:.2f}')
print('交互作用 ' + f'{QAB:.2f}'+' '*(15-len(f'{QAB:.2f}')) + f'{fab:d}'+' '*(12-len(f'{fab:d}'))
+ f'{QAB_bar:.2f}'+' '*(12-len(f'{QAB_bar:.2f}')) + f'{FAB:.2f}')
print('误差 ' + f'{QE:.2f}'+' '*(15-len(f'{QE:.2f}')) + f'{fe:d}'+' '*(12-len(f'{fe:d}'))
+ f'{QE_bar:.2f}'+' '*(12-len(f'{QE_bar:.2f}')))
print('总和 ' + f'{QT:.2f}'+' '*(15-len(f'{QT:.2f}')) + f'{ft:d}'+' '*(12-len(f'{ft:d}')))
print('=========================================================\n')
print("结论:")
'判断FA,FB是不是在拒绝域{FA≥F_alpha(fa,fe)},{FB≥F_alpha(fb,fe)},{FAB≥F_alpha(fab,fe)}里面'
FA_alpha = f.isf(alpha,fa,fe) # FA(fa,fe)的上α分位数
if FA>FA_alpha:
print(f"FA={FA:.2f} ≥ FA_{alpha:.2f}({fa:d}, {fe:d})={FA_alpha:.2f}，在显著性水平{alpha:.2f}上 拒绝 原假设H01。")
else:
print(f"FA={FA:.2f} < FA_{alpha:.2f}({fa:d}, {fe:d})={FA_alpha:.2f}，在显著性水平{alpha:.2f}上 接受 原假设H01。")
FB_alpha = f.isf(alpha,fb,fe) # FB(fb,fe)的上α分位数
if FB > FB_alpha:
print(f"FB={FB:.2f} ≥ FB_{alpha:.2f}({fb:d}, {fe:d})={FB_alpha:.2f}，在显著性水平{alpha:.2f}上 拒绝 原假设H02。")
else:
print(f"FB={FB:.2f} < FB_{alpha:.2f}({fb:d}, {fe:d})={FB_alpha:.2f}，在显著性水平{alpha:.2f}上 接受 原假设H02。")
FAB_alpha = f.isf(alpha,fab,fe) # FAB(fab,fe)的上α分位数
if FAB > FAB_alpha:
print(f"FAB={FAB:.2f} ≥ FAB_{alpha:.2f}({fab:d}, {fe:d})={FAB_alpha:.2f}，在显著性水平{alpha:.2f}上 拒绝 原假设H03。")
else:
print(f"FAB={FAB:.2f} < FAB_{alpha:.2f}({fab:d}, {fe:d})={FAB_alpha:.2f}，在显著性水平{alpha:.2f}上 接受 原假设H03。")
例2
if __name__=="__main__": #主函数
data3 = np.array([[[71,73],[72,73],[75,73],[77,75]],
[[73,75],[76,74],[78,77],[74,74]],
[[76,73],[79,77],[74,75],[74,73]],
[[75,73],[73,72],[70,71],[69,69]]])
alpha = 0.05
Two_way_ANOVA_repeated(data3, alpha)
作者：不想当程序员（☆∀☆）
