-
Notifications
You must be signed in to change notification settings - Fork 0
/
Multinomial.cpp
55 lines (48 loc) · 1.33 KB
/
Multinomial.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
#include <RcppArmadillo.h>
#include <omp.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using arma::uword;
arma::Col<int> rmultinom_1(const uword &size,
arma::vec &probs,
const uword &N) {
// ...
arma::Col<int> out(N);
rmultinom(size, probs.begin(), N, out.begin());
return out;
}
// [[Rcpp::export]]
arma::Mat<int> rmultinom_arma(const uword &n,
const uword &size,
arma::vec &probs) {
// ...
uword N = probs.size();
arma::Mat<int> sim(N, n);
for (uword i=0; i<n; ++i)
sim.col(i) = rmultinom_1(size, probs, N);
return sim;
}
// [[Rcpp::export]]
arma::Mat<int> rmultinom_omp(const uword &n,
const uword &size,
arma::vec &probs,
const int &cores=1) {
// ...
omp_set_num_threads(cores);
uword N = probs.size();
arma::Mat<int> sim(N, n);
#pragma omp parallel for
for (uword i=0; i<n; ++i)
sim.col(i) = rmultinom_1(size, probs, N);
return sim;
}
// [[Rcpp::export]]
IntegerVector rpois_arma(const uword &n,
const arma::vec &lambda) {
// ...
uword size = lambda.size();
IntegerVector sim(n);
for (uword i=0; i<n; ++i)
sim[i] = R::rpois(lambda[i%size]);
return sim;
}