用Scipy做ward聚类分析,但是他只输出聚类结果并未输出统计量,所以自己写了一个根据输出Z计算统计量的代码。Scipy库ward聚类算的distance和正常的离差平方和衡量的类间距离差了一个常数倍,需要处理一下,代码可能写的有点冗长,能用就行。
def Cluster_W(index,Z):
if index < Z.shape[0]+1:
W = 0
return W
else:
new_index = index - (Z.shape[0]+1)
k = int(Z[new_index,0])
l = int(Z[new_index,1])
W = Cluster_W(k,Z) + Cluster_W(l,Z) + Z[new_index,2]
return W
def Cluster_N(index,Z):
if index < Z.shape[0]+1:
N = 1
return N
else:
new_index = index - (Z.shape[0]+1)
N = int(Z[new_index,3])
return N
def WARD_Stat(Z):
## 根据Z计算每一步增加的离差平方和,输出统计量
d = (Z[:,2]**2)/4
T = d.sum()
W = np.insert(d,0,0)
P = 0
cluster = len(W)
Rk2 = np.zeros(len(W))
Partial_Rk2 = np.zeros(len(W))
pseudo_F = np.zeros(len(W))
pseudo_t = np.zeros(len(W))
N_Cluster = np.zeros(len(W))
for i in range(len(W)):
N_Cluster[i] = cluster
P += W[i]
Rk2[i] = 1-P/T
if i == 0:
pseudo_F[i] = np.nan
Partial_Rk2[i] = np.nan
pseudo_t[i] = np.nan
elif i == len(W)-1:
pseudo_F[i] = np.nan
Partial_Rk2[i] = Rk2[i-1]-Rk2[i]
k = int(Z[i-1,0])
l = int(Z[i-1,1])
Wk = Cluster_W(k,Z)
Wl = Cluster_W(l,Z)
pseudo_t[i] = ((Partial_Rk2[i]*T)*(Cluster_N(k,Z)+Cluster_N(l,Z)-2))/(Wk+Wl)
else:
Partial_Rk2[i] = Rk2[i-1]-Rk2[i]
pseudo_F[i] = ((T-P)/P)*((len(W)+1-cluster)/(cluster-1))
k = int(Z[i-1,0])
l = int(Z[i-1,1])
if k < len(W) and l < len(W):
pseudo_t[i] = np.nan
else:
Wk = Cluster_W(k,Z)
Wl = Cluster_W(l,Z)
pseudo_t[i] = ((Partial_Rk2[i]*T)*(Cluster_N(k,Z)+Cluster_N(l,Z)-2))/(Wk+Wl)
cluster -= 1
stat = pd.DataFrame(np.concatenate((N_Cluster.reshape(-1,1),Rk2.reshape(-1,1),Partial_Rk2.reshape(-1,1),pseudo_F.reshape(-1,1),pseudo_t.reshape(-1,1)),axis=1),columns=['N_Cluster','RK2','Partial_RK2','Pseudo_F','Pseudo_t2'])
return stat
最后输出一个dataframe,统计量相关的意义可自行查阅多元统计分析的书。