3

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

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

4

2 回答 2

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

干杯，戴维德

于 2013-09-23T19:26:00.327 回答
1

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

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

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

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

于 2013-03-20T02:32:48.690 回答