7

この質問は、さまざまなサンプル サイズと確率を持つ多項分布からの効率的なサンプリングに関するものです。以下に、私が使用したアプローチについて説明しますが、インテリジェントなベクトル化で改善できるかどうか疑問に思っています。

複数の集団間での生物の分散をシミュレートしています。集団からの個体は、確率 で集団に分散jします。母集団 1 の初期存在量が 10 で、母集団 1、2、3への分散確率がそれぞれ与えられると、次の式で分散プロセスをシミュレートできます。ip[i, j]c(0.1, 0.3, 0.6)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ここで、行 th 列の要素の値は、人口から人口へj移動する個体の数です。この行列の は、新しい人口サイズを示します。jirowSums

これを何度も繰り返したいのですが、確率行列は一定ですが、(事前定義された) 初期存在量は変化します。次の小さな例はこれを実現しますが、より大きな問題では非効率的です。結果のマトリックスは、集団の初期存在量が異なる5つのシミュレーションのそれぞれについて、3つの集団のそれぞれにおける分散後の存在量を示します。

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

この問題をより適切にベクトル化する方法はありますか?

4

3 に答える 3

7

これは多多項式の 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

gsl バージョンは、純粋な R バージョンよりも約 6 倍高速です。

于 2014-04-16T05:24:47.273 に答える