R中是否有Python bisect.bisect()的对应函数？

Question

我想从离散分布中抽样。

我有一个矩阵，叫做pi，里面包含了一些概率向量（每一行的列数相同，而且每一行的总和都是1）。

在Python中，我可以这样做：

cumsumpi = cumsum(pi, axis = 1)
[bisect.bisect(k, random.rand()) for k in cumsumpi]

这样就能根据pi给出的概率得到抽样的结果。

现在我想用R来实现这个功能。我知道R里面有一个叫“sample”的函数，但它似乎用的是不同的算法，所以即使我在两个地方都使用相同的set.seed()，得到的抽样结果也不一样。

我使用rpy2来确保在Python中得到和R中完全一样的随机抽样结果。例如，

我用的是：

[bisect.bisect(k, asarray(robjects.r('runif(1)'))) for k in cumsumpi]

请告诉我在R中有没有其他函数可以做同样的事情，而不是用sample。

-Joon

编辑：我用以下方法成功得到了完全一样的抽样结果，但速度比较慢。

    cumsumpi = t(apply(pi, 1, cumsum))

    getfirstindx = function(cumprobs) {
        return(which(cumprobs > runif(1))[1])
    }

    apply(cumsumpi, 1, getfirstindx)

Answer 1

我想找的是 findInterval - 查找区间的数字或索引。 :)

Answer 2

我看不太懂你问题的标题和内容之间的关系。不过，没关系，这里有一个R语言的函数，它和Python里的bisect功能是一样的：

这个gtools包里有一个二分查找的函数，叫做binsearch，它和Python的bisect几乎一模一样，比如：

# search for 25 in the range 0 through 100
> binseaerch(fun = function(x) x - 25, range=c(0, 100))

$call
binsearch(fun = function(x) x - 25, range = c(0, 100))

$numiter
[1] 2

$flag
[1] "Found"

$where
[1] 25

$value
[1] 0

Answer 3

这里有一种不同的方法，它不使用apply，而是通过向量化来完成操作。初步检查显示，这种方法的速度是之前的两倍，但还需要更详细的研究。

cumsumpi = t(apply(pi, 1, cumsum));
u = runif(nrow(cumsumpi));

max.col((cumsumpi > u) * 1, "first")

为了进一步加快速度，可以考虑对每一行计算累计列总和的操作进行向量化。如果你能运行一个性能分析工具来检查你的R代码，告诉我这个步骤是否是瓶颈。