当样本大小和概率变化时，高效的多项式采样

Question

这个问题涉及到在多项式分布中高效采样，其中样本大小和概率不同。以下是我所使用的方法，但我想知道是否可以通过一些智能矢量化来改进它。

我正在模拟生物在多个种群之间的扩散。来自种群j的个体以概率p[i,j]分散到种群i。假设在种群1中有10个个体，且分散概率为c(0.1,0.3,0.6)到分别到种群1、2、3，我们可以用rmultinom来模拟扩散过程：

set.seed(1)
rmultinom(1, 10, c(0.1, 0.3, 0.6))

#      [,1]
# [1,]    0
# [2,]    3
# [3,]    7

我们可以将其扩展到考虑 n 个源种群：

set.seed(1)
n <- 3
p <- replicate(n, diff(c(0, sort(runif(n-1)), 1)))
X <- sample(100, n)

上面，p是从一个种群（列）移动到另一个种群（行）的概率矩阵，X是初始种群大小的向量。现在可以使用以下方式模拟每对种群之间（以及留在原地的个体）之间的迁移数量：

sapply(seq_len(ncol(p)), function(i) {
  rmultinom(1, X[i], p[, i])  
})

#      [,1] [,2] [,3]
# [1,]   19   42   11
# [2,]    8   18   43
# [3,]   68    6    8

在第i行和第j列的元素值是从种群j移动到种群i的个体数量。该矩阵的rowSums给出了新的种群大小。

我希望用不同（预定义的）初始丰度多次重复此过程，但概率矩阵保持不变。以下小例子可以实现此目标，但在处理更大问题时效率较低。所得到的矩阵给出了拥有不同初始丰度的五个种群中每个种群的扩散后丰度。

X <- matrix(sample(100, n*5, replace=TRUE), nrow=n)

apply(sapply(apply(X, 2, function(x) {
  lapply(seq_len(ncol(p)), function(i) {
    rmultinom(1, x[i], p[, i])  
  })
}), function(x) do.call(cbind, x), simplify='array'), 3, rowSums)

#      [,1] [,2] [,3] [,4] [,5]
# [1,]   79   67   45   28   74
# [2,]   92   99   40   19   52
# [3,]   51   45   16   21   35

有没有更好的向量化这个问题的方法？

- jbaums

---

照常，快速回答是：你尝试过使用“Rcpp”吗？请尝试访问http://www.gnu.org/software/gsl/manual/html_node/The-Multinomial-Distribution.html。 - Gary Weissman

@GaryWeissman：还没有——我很好奇在冒险之前是否可能进行纯R改进。 - jbaums

谢谢链接，@GaryWeissman - 你有没有想法，是否使用Rcpp中的Environment stats("package:stats"); Function rmultinom = stats["rmultinom"];比在这里（http://lists.r-forge.r-project.org/pipermail/rcpp-devel/2011-June/002393.html）所示的方式更快/更好？ - jbaums

我刚刚测试了你提到的示例，并将其与本地rmultinom函数进行比较，后者快了一个数量级。除非Rcpp可以比您的嵌套应用程序更快地循环遍历不同的总体，否则加速可能并不那么明显。 - Gary Weissman

我了解您已经有足够的计算答案了。但是请记住，如果您遇到速度瓶颈，最简单的方法可能是采用分析路线。您所描述的是离散马尔可夫链，因此可以快速获得平衡概率。此外，对于大N，多项式分布几乎是多元正态分布，可以通过这种方式进行近似。因此，如果速度较慢，我会选择算法实现而不是算法实现的编码，因为这样您的问题可能更容易回答。 - evolvedmicrobe

Answer 1

这是一个RcppGSL实现的多重多项式分布。然而，它需要您单独安装gsl...这可能不是非常实际。

// [[Rcpp::depends(RcppGSL)]]

#include <RcppGSL.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <unistd.h>            // getpid

Rcpp::IntegerVector rmn(unsigned int N, Rcpp::NumericVector p, gsl_rng* r){

    size_t K = p.size();

    Rcpp::IntegerVector x(K);
    gsl_ran_multinomial(r, K, N, p.begin(), (unsigned int *) x.begin());
    return x;             // return results vector
}

Rcpp::IntegerVector gsl_mmm_1(Rcpp::IntegerVector N, Rcpp::NumericMatrix P, gsl_rng* r){
    size_t K = N.size();
    int i;
    Rcpp::IntegerVector x(K);
    for(i=0; i<K; i++){
        x += rmn(N[i], P(Rcpp::_, i), r);
    }
    return x;
}

// [[Rcpp::export]]
Rcpp::IntegerMatrix gsl_mmm(Rcpp::IntegerMatrix X_, Rcpp::NumericMatrix P){
    int j;
    gsl_rng * r = gsl_rng_alloc (gsl_rng_mt19937);
    long seed = rand()/(((double)RAND_MAX + 1)/10000000) * getpid();
    gsl_rng_set (r, seed);
    Rcpp::IntegerMatrix X(X_.nrow(), X_.ncol());
    for(j=0; j<X.ncol(); j++){
        X(Rcpp::_, j) = gsl_mmm_1(X_(Rcpp::_,j), P, r);
    }
    gsl_rng_free (r);
    return X;
}

我还将其与纯R实现和jbaums版本进行了比较

library(Rcpp)
library(microbenchmark)
sourceCpp("gsl.cpp")

P = matrix(c(c(0.1,0.2,0.7),c(0.3,0.3,0.4),c(0.5,0.3,0.2)),nc=3)
X = matrix(c(c(30,40,30),c(20,40,40)), nc=2)

mmm = function(X, P){
    n = ncol(X)
    p = nrow(X)
    Reduce("+", lapply(1:p, function(j) {
        Y = matrix(0,p,n)
        for(i in 1:n) Y[,i] = rmultinom(1, X[j,i], P[,j])
        Y
    }))
}

jbaums = function(X,P){
    apply(sapply(apply(X, 2, function(x) {
      lapply(seq_len(ncol(P)), function(i) {
        rmultinom(1, x[i], P[, i])
      })
    }), function(x) do.call(cbind, x), simplify='array'), nrow(X), rowSums)
}
microbenchmark(jbaums(X,P), mmm(X,P), gsl_mmm(X, P))

这是结果

> microbenchmark(jbaums(X,P), mmm(X,P), gsl_mmm(X, P))
Unit: microseconds
          expr     min       lq  median       uq     max neval
  jbaums(X, P) 165.832 172.8420 179.185 187.2810 339.280   100
     mmm(X, P)  60.071  63.5955  67.437  71.5775  92.963   100
 gsl_mmm(X, P)  10.529  11.8800  13.671  14.6220  40.857   100

比我的纯R版本快了约6倍。

- Randy Lai

---

我在 gsl 版本中发现了一个 bug，现在已经修复了。 - Randy Lai

头文件的路径是否只是相对于 gsl.cpp 所在的目录？它们只是在文件顶部提供，还是我需要做其他事情来注册它们到 Rcpp 中？这里是相对新手。 - jbaums

不，它们应该在系统PATH变量列出的某个地方。当您安装gsl时，相应的位置应该被添加到PATH中。如果您可以安装RcppGSL，那么意味着Rcpp可以找到gsl。 - Randy Lai

我还没有能够在Windows上安装gsl并测试它，但是你的基准测试表明它将是合适的。 - jbaums

如果您无法使用gsl，请考虑使用我的R实现。 - Randy Lai

谢谢，是的，“Reduce”方法也很有用。看看每个方法在多项试验数量增加时如何扩展会很有趣。我下面的C++版本是对您的“Reduce”方法的一点改进，但是随着试验次数的增加，它很快就比我最初的R版本变慢了。 - jbaums

Answer 2

例如：

# make the example in Rcpp you mention:
library(Rcpp)
library(inline)
src <- 'Environment stats("package:stats");
Function rmultinom = stats["rmultinom"];
NumericVector some_p(1000, 1.0/1000);
return(rmultinom(1,1, some_p));'

fx <- rcpp(signature(), body=src)

# now compare the two
library(rbenchmark)
benchmark(fx(),rmultinom(1,1,c(1000,1/1000)),replications=10000)

#                            test replications elapsed relative user.self sys.self user.child sys.child
#    1                       fx()        10000   1.126   13.901     1.128        0          0         0
#    2 rmultinom(1, 1, c(1/1000))        10000   0.081    1.000     0.080        0          0         0

- Gary Weissman

---

在C环境下运行R函数通常比在R环境下运行R函数慢得多。 - Randy Lai

感谢这个基准测试，也感谢您的提示，@RandyLai。 Gary - 您指出的这个示例如何使用Rcpp实现？ - jbaums

我似乎无法将GSL库链接起来，但代码看起来应该是这样的：fx2 <- cxxfunction(signature(), body='return(gsl_ran_multinomial(1,1,1000,1.0/1000))', includes='#include <gsl/gsl_ran_dist.h>', plugin='RcppGSL') - Gary Weissman

@jbaums，你知道如何安装gsl吗？因为RcppGSL不带gsl库，所以你需要自己安装。 - Randy Lai

@RandyLai：谢谢提醒。我还没有尝试过...在Windows上似乎是可能的，但听起来需要通过Cygwin运行？如果是这样，我想R终端（或更好的是RStudio）也必须通过Cygwin运行...这样会行吗？虽然理想情况下我正在寻找一个简单的跨平台解决方案。 - jbaums

Answer 3

我发现BH包将boost库引入了该项目中。这使得以下操作成为可能，其产生的输出与@RandyLai的gsl_mmm以及我上面提出的代码相同。（我相信启用c++11支持应该可以在没有BH的情况下使用random。）

// [[Rcpp::depends(BH)]]
#include <Rcpp.h>

#include <boost/random.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/discrete_distribution.hpp>

using namespace Rcpp;

typedef boost::mt19937 RNGType;
RNGType rng(123);


NumericVector rowSumsC(IntegerMatrix x) {
  int nrow = x.nrow(), ncol = x.ncol();
  IntegerVector out(nrow);

  for (int i = 0; i < nrow; i++) {
    double total = 0;
    for (int j = 0; j < ncol; j++) {
      total += x(i, j);
    }
    out[i] = total;
  }
  return wrap(out);
}

// [[Rcpp::export]]
IntegerMatrix rmm(IntegerMatrix X, NumericMatrix P) {
  int niter = X.ncol(), nx = X.nrow();
  IntegerMatrix out(nx, niter);
  for (int j = 0; j < niter; j++) {
    IntegerMatrix tmp(nx, nx);
    for (int i = 0; i < nx; i++) {
      for (int n = 0; n < X(i, j); n++) {
        boost::random::discrete_distribution<> dist(P(_, i));
        tmp(dist(rng), i)++;
      }
    }
    out(_, j) = rowSumsC(tmp);
  }
  return out;
}

rowSumsC由@hadley提供，此处。

然而，在我的机器上，这比Randy的gsl_mmm慢得多，而且在有许多试验时比我的R版本还要慢。我怀疑这是由于编码效率低下，但是boost的discrete_distribution也会单独执行每个多项式试验，而使用gsl时，该过程似乎是向量化的。我对c++很陌生，不确定是否可以更加高效。

P <- matrix(c(c(0.1, 0.2, 0.7), c(0.3, 0.3, 0.4), c(0.5, 0.3, 0.2)), nc=3)
X <- matrix(c(c(30, 40, 30), c(20, 40, 40)), nc=2)
library(BH)
microbenchmark(jbaums(X, P), rmm(X, P))

# Unit: microseconds
#          expr     min       lq  median       uq     max neval
#  jbaums(X, P) 124.988 129.5065 131.464 133.8735 348.763   100
#     rmm(X, P)  59.031  60.0850  62.043  62.6450 117.459   100

- jbaums