python WARD聚类并输出统计量

用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，统计量相关的意义可自行查阅多元统计分析的书。