Python和C++及MATLAB距离相关性生物医学样本统计量算法及数据科学

🎯要点

统计观测值之间距离计算
代谢组学和脂质组学分析
相关距离矩阵计算
卡方检验偏差校正
快速计算距离协方差算法
大规模生物系统分析
距离矩阵相关性测试
石油勘探统计学关系

Python距离矩阵

在数学、计算机科学，尤其是图论中，距离矩阵是一个方阵（二维数组），其中包含一组元素之间成对的距离。根据所涉及的应用，用于定义此矩阵的距离可能是也可能不是度量。如果有 N 个元素，则此矩阵的大小为 N×N。在图论应用中，元素通常被称为点、节点或顶点。

一般来说，距离矩阵是某个图的加权邻接矩阵。在网络（即为弧分配权重的有向图）中，网络两个节点之间的距离可以定义为连接两个节点的最短路径上权重之和的最小值（其中路径中的步数是有界的）。[2] 这个距离函数虽然定义明确，但不是度量。除了需要能够组合和比较权重之外，不需要对权重进行任何限制，因此在某些应用中会使用负权重。由于路径是有向的，因此无法保证对称性，并且如果存在负权重循环，距离矩阵可能不是空心的（并且如果没有对步数的限制，矩阵可能未定义）。

上述的代数公式可以通过使用最小加代数来获得。该系统中的矩阵乘法定义如下：给定两个 $\times n$ 矩阵 $A=\left(a_{i j}\right)$ 和 $B=\left(b_{i j}\right)$ ，它们距离积 $C=\left(c_{i j}\right)=A \star B$ 定义为 $\times n$ 矩阵，使得
$c_{i j}=\min _{k=1}^n\left\{a_{i k}+b_{k j}\right\}$
请注意，非直接连接的对角线元素需要设置为无穷大或合适的大值，以使 min-plus 操作正常工作。这些位置的零将被错误地解释为没有距离、成本等的边。

Python示例说明如何对距离矩阵进行排序，以便最终（分层）集群是明显的。

import numpy as np
from scipy.spatial.distance import pdist, squareform
from sklearn import datasets
import matplotlib.pyplot as plt

iris = datasets.load_iris()
iris.data.shape

(150, 4)

该数据集包含 R^4 中的 150 个点（由 4 个特征描述的 150 朵花）。使用 pdist 我们可以计算一个 150 x 150 的距离矩阵，如下所示。

dist_mat = squareform(pdist(iris.data))

N = len(iris.data)
plt.pcolormesh(dist_mat)
plt.colorbar()
plt.xlim([0,N])
plt.ylim([0,N])
plt.show()

我们可以看到有 3 个明显的簇：一个相当密集（左下）并且距离其他两个簇很远，另外两个非常接近，但与第三个簇（左下）的距离不同。

这些聚类对应于可以在 iris.target 中找到的 3 个类：一个是线性可分的（下面的 0 类，对应于距离矩阵中左下方的聚类），另外两个不是线性可分的（类 1 和 2，对应于最右边的两个聚类）。

iris.target

array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2])

为了展示所提出的相异矩阵排序方法的有用性，我们对数据点应用随机排列。这是相关的，因为数据点通常以任意顺序提供并且通常没有任何标签。

X = iris.data[np.random.permutation(N),:]

dist_mat = squareform(pdist(X))

plt.pcolormesh(dist_mat)
plt.xlim([0,N])
plt.ylim([0,N])
plt.show()

在该过程结束时，所有点都根据树状图施加的顺序进行排序。

def seriation(Z,N,cur_index):

    if cur_index < N:
        return [cur_index]
    else:
        left = int(Z[cur_index-N,0])
        right = int(Z[cur_index-N,1])
        return (seriation(Z,N,left) + seriation(Z,N,right))
    
def compute_serial_matrix(dist_mat,method="ward"):

    N = len(dist_mat)
    flat_dist_mat = squareform(dist_mat)
    res_linkage = linkage(flat_dist_mat, method=method,preserve_input=True)
    res_order = seriation(res_linkage, N, N + N-2)
    seriated_dist = np.zeros((N,N))
    a,b = np.triu_indices(N,k=1)
    seriated_dist[a,b] = dist_mat[ [res_order[i] for i in a], [res_order[j] for j in b]]
    seriated_dist[b,a] = seriated_dist[a,b]
    
    return seriated_dist, res_order, res_linkage