协方差矩阵的齐性检验（接上一篇文章）

关于协方差同质性检验，我也是一知半解，不过多讲解，自己也很懵。
本文讲述对于两分类问题的协方差矩阵检验，和多分类的协方差矩阵的检验

两分类的协方差矩阵齐性检验

Σ1是类别1的协方差，Σ2是类别2的协方差，Σ是两个协方差的联合协方差（即图中的Σ_hat，和S）。
式子中 tr 表示trace，即沿着对角线求和。
p是维数，即数据有几个特征

因为在贝叶斯判别的式子中如果两分类协方差相等，那么用的是他们的联合协方差做计算。

所以这个检验方法是检验Σ1 和联合协方差的差异性，和Σ2与联合协方差的差异性。
如果他们的差异性都不显著，则说明他们的协方差齐性，可以用联合协方差计算。

还是 用这个例子的数据。

x = np.array([[1.14,1.78],[1.18,1.96],[1.20,1.86],[1.26,2.00],[1.28,2.00],[1.30,1.96],[1.24,1.72],[1.36,1.74],[1.38,1.64],
              [1.38,1.82],[1.38,1.90],[1.40,1.70],[1.48,1.82],[1.54,1.82],[1.56,2.08]])

y = np.array([1,1,1,1,1,1,2,2,2,2,2,2,2,2,2])

import numpy as np
import scipy.stats as stats

# 把数据分类分开
x1 = x[y==1]
x2 = x[y==2]

def cov_homogeneity_2(arr1,arr2):
	'''
	两个类别的协方差齐性检验，这里返回True和False，以及Q1,Q2，和拒绝阈
	arr1 ，arr2分别是两个类别的数据特征
	'''
	# p 是维数，数据有几维，也就是几个特征
	if (arr1.shape[1] == arr2.shape[1]):
        p = arr1.shape[1]
    else:raise(ValueError("数据维度不一样"))
 
    a1 = np.array(arr1,dtype=np.float64)
    a2 = np.array(arr2,dtype=np.float64)
    
    # 计算样本个数，和协方差，以及联合协方差
    n1,n2 = a1.shape[0],a2.shape[0]
    s1 = np.cov(a1,rowvar=False)
    s2 = np.cov(a2,rowvar=False)
    s = ((n1-1)*s1+(n2-1)*s2)/(n1+n2-2)
    
    
    
    # 按照公式计算Q的统计值
    Q1 = (n1-1)*(np.log(np.linalg.det(s)) - np.log(np.linalg.det(s1)) - p + 
    	np.matrix.trace(np.linalg.inv(s).dot(s1)))
    	
    Q2 = (n2-1)*(np.log(np.linalg.det(s)) - np.log(np.linalg.det(s2)) - p + 
    	np.matrix.trace(np.linalg.inv(s).dot(s2)))
    df = p*(p+1)/2
    # 因为这是卡方分布，先计算临界值，也就是拒绝阈
    lingjie = stats.chi2.ppf(0.95,df)

	# 如果都比拒绝阈小，说明两个协方差无显著差异（H0：协方差无显著差异）
    if (Q1<lingjie) and (Q2 < lingjie):
        return True,Q1,Q2,lingjie
    else:return False,Q1,Q2,lingjie

print(cov_homogeneity_2(x1,x2))
(True, 2.5784296441237355, 0.7417513518834973, 7.814727903251179)
# 返回True，所以这个数据的协方差通过同质性检验，可以认为协方差相等

多分类的协方差矩阵齐性检验

先讨论三个总体的情况，其实三个总体也可以用上面的方法两两比较，比较多次就可得出结论。同样的，四个五个也是一样，多比较几次总会有结果。
下面看图：

上图是三个类别，如果把他推广到k个类别，唯一不同的，就是联合协方差矩阵S
如果是k个类别，那么

其中 n1 一直加到 nk 就是总样本 n 个。
这里注意当n=k时，比如三个样本三个类别，那么这时候也无需分类，因为数据太少，判别不了。。所以每个类别至少有两个数据，才能进行判别。即n>k
直接给出推广到k个类别的判断函数

def cov_homogeneity_n(x,y):
    if x.shape[0]!=y.shape[0]:
    	return "x和y数据不是相同的个数"
    
    n = x.shape[0] # 样本的总个数
    p = x.shape[1] # 维度
    k = np.unique(y).size # 样本的类别数
    if n=k:return "请检查数据"
    
    s1=0 # 联合协方差矩阵用的
    s2=0 # 这表示M后面减去部分的求和项，也在循环中求出
    df=p*(p+1)*(k-1)/2 # 自由度
    
    # 把d分成前面固定的数和后面求和项两部分的乘积，在循环内部算出求和项d2
    d1 = (2*p**2+3*p-1)/(6*(p+1)*(k-1))
    d2=0
    
    for y_i in np.unique(y):
        x_i = x[y==y_i]
        n_y_i = x_i.shape[0] # 属于这一类的样本个数
        
        d2+=(1/(n_y_i-1)-1/(n-k))
        
        s_i = np.cov(x_i,rowvar=False)
        
        s1+=(n_y_i-1)*s_i/(n-k) # 联合协方差矩阵，其实也是求和项，分母相同，分子是每一项的求和
        
        s2+=(n_y_i-1)*np.log(np.linalg.det(s_i))
    
    # 最终的公式：
    d = d1*d2
    
    S = s1
    
    M = (n-k)*np.log(np.linalg.det(S))-s2
    
    T = (1-d)*M
    
    lingjie = stats.chi2.ppf(0.95,df)
    if T<lingjie:
        return True,T,lingjie
    else:return False,T,lingjie

给出下列数据

x = np.array([[261.01,7.36],[185.39,5.99],[249.58,6.11],[137.13,4.35],[231.34,8.79],[231.38,8.53],[260.25,10.02],
 [259.51,9.79],[273.84,8.79],[303.59,8.53],[231.03,6.15],[308.90,8.49],[258.69,7.16],[355.54,9.43],
 [476.69,11.32],[316.12,8.17],[274.57,9.67],[409.42,10.49],[330.34,9.61],[331.47,13.72],[352.50,11.00],
 [347.31,11.19],[189.59,5.46]])

y = np.array([1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,3,3,3,3,3])

cov_homogeneity_n(x,y)

out:
(True, 8.138112152750503, 12.591587243743977)

这组数据三个类别中协方差矩阵齐性，即无显著性差异。

SPSS中判断协方差的方法

在spss中也提供了判断协方差齐性的假设检验。
还是上面的数据，三个类别，spss中在分析-分类-判别式，点进去就是spss的判别分析，这里只针对说明协方差的假设检验。

在上图统计按钮点开后得到下面：

box’M检验，就是协方差的假设性检验，H0为协方差相等。

这里用到的是F检验，具体的检验方法我没搞明白，网上资料有限，知道的可以留言给我，谢谢了。

结果也是这三个类别的协方差相等，没有理由拒绝原假设。

---

另外，顺嘴一提，上图统计的对话框里，函数系数列表中有Fisher和未标准化两种选择，

其中注意Fisher选项，为贝叶斯判别，是用概率来计算属于哪一个类别
他给出的系数：

有几个类别就会有几组函数的系数，
要判断一个数据时，把数据分别带入三个方程，计算得到的结果为y1,y2,y3
那么属于第一类的概率为

以此计算属于其他类的概率，哪个概率大就判给谁

而未标准化则给出Fisher判别的未标准化系数，因为默认给出标准化的
由Fisher判别公式计算离质心的距离来判断属于哪一类别。
这里容易混淆，所以记一下！