在 Numpy/Python 中快速稀疏矩阵答案

【问题标题】：Quickly rarefy a matrix in Numpy/Python在 Numpy/Python 中快速稀疏矩阵
【发布时间】：2013-03-19 18:47:29
【问题描述】：

我需要（快速）稀疏矩阵。

Rarefaction - 将丰度矩阵转换为均匀的采样深度。

在此示例中，每一行都是一个样本，采样深度是该行的总和。我想通过min(rowsums(matrix)) 样本随机采样（替换）矩阵。

假设我有一个矩阵：

>>> m = [ [0, 9, 0],
...       [0, 3, 3],
...       [0, 4, 4] ]

稀疏函数逐行随机抽样，替换为min(rowsums(matrix)) 次（在本例中为 6）。

>>> rf = rarefaction(m)
>>> rf
    [ [0, 6, 0],  # sum = 6
      [0, 3, 3],  # sum = 6
      [0, 3, 3] ] # sum = 6

结果是随机的，但行总和始终相同。

>>> rf = rarefaction(m)
>>> rf
    [ [0, 6, 0],   # sum = 6
      [0, 2, 4],   # sum = 6
      [0, 4, 2], ] # sum = 6

PyCogent 有一个函数可以逐行执行此操作，但在大型矩阵上速度非常慢。

我感觉 Numpy 中有一个函数可以做到这一点，但我不确定它会被调用什么。

【问题讨论】：

我认为nowsums 真的意味着rowsums？

标签： python numpy

【解决方案1】：

import numpy as np
from numpy.random import RandomState

def rarefaction(M, seed=0):
    prng = RandomState(seed) # reproducible results
    noccur = np.sum(M, axis=1) # number of occurrences for each sample
    nvar = M.shape[1] # number of variables
    depth = np.min(noccur) # sampling depth

    Mrarefied = np.empty_like(M)
    for i in range(M.shape[0]): # for each sample
        p = M[i] / float(noccur[i]) # relative frequency / probability
        choice = prng.choice(nvar, depth, p=p)
        Mrarefied[i] = np.bincount(choice, minlength=nvar)

    return Mrarefied

例子：

>>> M = np.array([[0, 9, 0], [0, 3, 3], [0, 4, 4]])
>>> M
array([[0, 9, 0],
       [0, 3, 3],
       [0, 4, 4]])
>>> rarefaction(M)
array([[0, 6, 0],
       [0, 2, 4],
       [0, 3, 3]])
>>> rarefaction(M, seed=1)
array([[0, 6, 0],
       [0, 4, 2],
       [0, 3, 3]])
>>> rarefaction(M, seed=2)
array([[0, 6, 0],
       [0, 3, 3],
       [0, 3, 3]])

干杯，大卫

【讨论】：

【解决方案2】：

我认为这个问题并不完全清楚。我想稀疏矩阵为您提供了您从原始矩阵的每个系数中获取的样本数？

查看链接中的代码，可能会加快速度。对转置矩阵进行操作并重写链接代码以对列而不是行进行操作。因为这将允许您的处理器更好地缓存它采样的值，即内存中的跳转更少。

其余的就像我会做的那样，使用 numpy （不一定意味着这是最有效的方法）。

如果您需要更快的速度，您可以尝试用 C++ 编写函数并使用 scipy.weave 将其包含到您的 python 中。在 C++ 中，我将查找每一行并构建一个 >0 位置的查找表，在等于查找表中项目数的范围内生成 min(rowsums(matrix)) 整数。我会累积查找表中每个位置的绘制频率，然后将这些数字放回数组中的正确位置。该代码应该只有几行代码。

【讨论】：