在 Numpy/Python 中快速稀疏矩阵

Question

我需要（快速）稀薄一个矩阵。

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 中有一个函数可以做到这一点，但我不确定它会被称为什么。

Answer 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]])

干杯，大卫