numpy 相当于 matlab dummyvar

Question

为了很好地处理类别变量，matlab 的 dummyvar 函数最pythonic的等价物是什么？

这是一个说明我的问题的示例，它使用 NxM 矩阵表示将 N 个数据点划分为

>> partitions
array([[1, 1, 2, 2, 1, 2, 2, 2, 1, 1],
   [1, 2, 2, 1, 2, 1, 2, 2, 2, 1],
   [1, 1, 1, 2, 2, 2, 1, 3, 3, 2]])

任务是高效统计任意两个数据点被归为同一类别的次数，并将结果存储在一个 NxN 矩阵中。在 matlab 中，这可以通过 dummyvar 的单行来实现，它为每个分区的每个类别创建一个列变量。

>> dummyvar(partitions)*dummyvar(partitions)'
ans =
 3     2     1     1     1     1     1     0     1     2
 2     3     2     0     2     0     2     1     2     1
 1     2     3     1     1     1     3     2     1     0
 1     0     1     3     1     3     1     1     0     2
 1     2     1     1     3     1     1     1     2     2
 1     0     1     3     1     3     1     1     0     2
 1     2     3     1     1     1     3     2     1     0
 0     1     2     1     1     1     2     3     2     0
 1     2     1     0     2     0     1     2     3     1
 2     1     0     2     2     2     0     0     1     3

我能想到的解决此任务的最有效方法是编写一个 O(n*m) 循环来模拟 dummyvar 的行为。 （请注意，下面的代码更喜欢partition.shape[0] partition.shape[1]，这通常是正确的，但假设是不安全的）。

dv=np.zeros((0,10))
for row in partitions:
  for val in xrange(1,np.max(row)+1):
    dv=np.vstack((dv,row==val))
np.dot(dv.T,dv)

当然，因为循环中的 vstack 效率非常低，这可以通过找到所需的大小并创建数组来改进，但我真的在寻找一个像在 matlab 中一样的单行器。

编辑：有关我正在做什么的更多信息，只是为了添加上下文。我正在用python（不存在python实现）编写用于分析大脑网络的库的库函数。现有的工作 matlab 源是可用的。由于特定领域的限制，输入的最大大小大约是几千个节点的网络。但是，基本上我编写的所有函数都必须很好地扩展到大型输入。

Answer 1

你可以做一点广播魔法来快速获得你的虚拟数组：

>>> partitions = np.array([[1, 1, 2, 2, 1, 2, 2, 2, 1, 1],
...                        [1, 2, 2, 1, 2, 1, 2, 2, 2, 1],
...                        [1, 1, 1, 2, 2, 2, 1, 3, 3, 2]])
>>> n = np.max(partitions)
>>> d = (partitions.T[:, None, :] == np.arange(1, n+1)[:, None]).astype(np.int)
>>> d = d.reshape(partitions.shape[1], -1)
>>> d.dot(d.T)
array([[3, 2, 1, 1, 1, 1, 1, 0, 1, 2],
       [2, 3, 2, 0, 2, 0, 2, 1, 2, 1],
       [1, 2, 3, 1, 1, 1, 3, 2, 1, 0],
       [1, 0, 1, 3, 1, 3, 1, 1, 0, 2],
       [1, 2, 1, 1, 3, 1, 1, 1, 2, 2],
       [1, 0, 1, 3, 1, 3, 1, 1, 0, 2],
       [1, 2, 3, 1, 1, 1, 3, 2, 1, 0],
       [0, 1, 2, 1, 1, 1, 2, 3, 2, 0],
       [1, 2, 1, 0, 2, 0, 1, 2, 3, 1],
       [2, 1, 0, 2, 2, 2, 0, 0, 1, 3]])

有一个明显的缺点是，即使一行只有几个不同的值，我们创建的虚拟数组也会为该行拥有尽可能多的列，就像值最多的行一样。但除非你有巨大的数组，否则它可能会比任何其他方法都快。

---

好吧，如果您想要一个可扩展的解决方案，您想为您的虚拟矩阵使用一个稀疏数组。如果您不熟悉 CSR 稀疏格式的细节，以下代码可能难以理解：

import scipy.sparse as sps
def sparse_dummyvar(partitions):
    num_rows = np.sum(np.max(partitions, axis=1))
    nnz = np.prod(partitions.shape)
    as_part = np.argsort(partitions, axis=1)
    # You could get s_part from the indices in as_part, left as
    # an exercise for the reader...
    s_part = np.sort(partitions, axis=1)
    mask = np.hstack(([[True]]*len(items_per_row),
                      s_part[:, :-1] != s_part[:, 1:]))
    indptr = np.where(mask.ravel())[0]
    indptr = np.append(indptr, nnz)

    return sps.csr_matrix((np.repeat([1], nnz), as_part.ravel(), indptr),
                          shape=(num_rows, partitions.shape[1],))

这将返回dummyvar(partitions) 的转置。您可以通过调用csc_matrix 而不是csr_matrix 并交换形状值来获得数组而无需转置。但是由于您只在矩阵及其转置的乘积之后，并且 scipy 在相乘之前将所有内容转换为 CSR 格式，因此它可能会稍微快一些。你现在可以这样做：

>>> dT = sparse_dummyvar(partitions)
>>> dT.T.dot(dT)
<10x10 sparse matrix of type '<type 'numpy.int32'>'
    with 84 stored elements in Compressed Sparse Column format>
>>> dT.T.dot(dT).A
array([[3, 2, 1, 1, 1, 1, 1, 0, 1, 2],
       [2, 3, 2, 0, 2, 0, 2, 1, 2, 1],
       [1, 2, 3, 1, 1, 1, 3, 2, 1, 0],
       [1, 0, 1, 3, 1, 3, 1, 1, 0, 2],
       [1, 2, 1, 1, 3, 1, 1, 1, 2, 2],
       [1, 0, 1, 3, 1, 3, 1, 1, 0, 2],
       [1, 2, 3, 1, 1, 1, 3, 2, 1, 0],
       [0, 1, 2, 1, 1, 1, 2, 3, 2, 0],
       [1, 2, 1, 0, 2, 0, 1, 2, 3, 1],
       [2, 1, 0, 2, 2, 2, 0, 0, 1, 3]])