Quickly rarefy a matrix in Numpy/Python

numpy, python

Solution

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

Example:

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

Cheers, Davide

Problem

I need to (quickly) rarefy a matrix. Rarefaction - transform abundance matrices to even sampling depth. In this example, each row is a sample and the sampling depth is the sum of the row. I want to randomly sample (with replacement) the matrix by `min(rowsums(matrix))` samples. Suppose I have a matrix: ``` >>> m = [ [0, 9, 0], ... [0, 3, 3], ... [0, 4, 4] ] ``` The rarefaction function goes row by row randomly sampling with replacement `min(rowsums(matrix))` times (which is 6 in this case). ``` >>> rf = rarefaction(m) >>> rf [ [0, 6, 0], # sum = 6 [0, 3, 3], # sum = 6 [0, 3, 3] ] # sum = 6 ``` The results are random but the row sums are always the same. ``` >>> rf = rarefaction(m) >>> rf [ [0, 6, 0], # sum = 6 [0, 2, 4], # sum = 6 [0, 4, 2], ] # sum = 6 ``` PyCogent has a function that does this row by row however it is very slow on large matrices. I have a feeling that there is a function in Numpy that can do this but I'm not sure what it would be called.

Original source