我正在将一个函数从R转换为Rcpp,并且已经苦苦挣扎了一段时间。我看过这样的一些材料一个,但我在RCPP一个初学者,我一直没能找出我目前做错了。
我在Rcpp中的最佳成绩如下,
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
Rcpp::NumericVector aux_pred_C(IntegerVector m,
NumericVector a,
NumericVector b,
int n_group,
IntegerVector group_index,
NumericMatrix theta_sample,
IntegerMatrix prior_parm,
int n_mcmc){
NumericVector p(n_mcmc);
group_index = group_index - 1;
for (int j = 0; j < n_mcmc; ++j) {
NumericMatrix parm(2, n_group);
NumericMatrix sample(n_mcmc, n_group);
for(IntegerVector::iterator i = group_index.begin(); i != group_index.end(); ++i) {
NumericVector temp = Rcpp::rbinom(m(i) , 1, theta_sample(j, i));
parm(0,i) = prior_parm(0,i) + a(i) + sum(temp);
parm(1,i) = prior_parm(1,i) + b(i) + m(i) - sum(temp);
sample(_,i) = Rcpp::rbeta(n_mcmc, parm(0,i), parm(1,i));
}
int i1 = group_index[0];
int i2 = group_index[1];
LogicalVector V = sample(_,i2) > sample(_,i1);
IntegerVector res = ifelse(V, 1, 0);
p(j) = mean(res);
}
return(p);
}
在第22、24、25和26行重复我的错误消息:
sampler_predictive_distribution.cpp:22:44: error: no match for call to '(Rcpp::IntegerVector {aka Rcpp::Vector<13>}) (Rcpp::traits::storage_type<13>::type*&)'
invalid conversion from 'Rcpp::Vector<13>::iterator' {aka 'int*'} to 'size_t' {aka 'long long unsigned int'} [-fpermissive] NumericVector temp = Rcpp::rbinom(m(i) , 1, theta_sample(j, i));
我的猜测是我没有正确使用迭代器。我遵循了这个例子。
最后,我在R中的可复制示例,
n.group <- 4
group.index <- c(1, 3)
m <- rep(NA, n.group)
m[group.index[1]] <- 50
m[group.index[2]] <- 50
a <- rep(NA, n.group)
a[group.index[1]] <- 20
a[group.index[2]] <- 20
b <- rep(NA, n.group)
b[group.index[1]] <- 30
b[group.index[2]] <- 30
n.mcmc <- 100
theta.sample.e1 <- matrix(NA, nrow = n.mcmc, ncol = n.group)
theta.sample.e1[, group.index[1]] <- rbeta(n.mcmc, a[group.index[2]], b[group.index[1]])
theta.sample.e1[, group.index[2]] <- rbeta(n.mcmc, a[group.index[2]], b[group.index[2]])
prior.parm.e1 <- matrix(1, ncol = 4, nrow = 2)
aux_pred <- function(m, a, b, n.group, group.index, theta.sample, prior.parm, n.mcmc){
p <- rep(NA, n.mcmc)
for (j in 1:n.mcmc){
parm <- matrix(NA, ncol = n.group, nrow = 2)
sample <- matrix(NA, ncol = n.group, nrow = n.mcmc)
for (i in group.index){
temp <- rbinom(m[i], size = 1, prob = theta.sample[j, i])
parm[1, i] <- prior.parm[1, i] + a[i] + sum(temp)
parm[2, i] <- prior.parm[2, i] + b[i] + length(temp) - sum(temp)
sample[, i] <- rbeta(n.mcmc, parm[1, i], parm[2, i])
}
p[j] <- mean(sample[, group.index[2]] - sample[, group.index[1]] > 0)
}
return(p)
}
aux_pred(m, a, b, n.group, group.index, theta.sample.e1, prior.parm.e1, n.mcmc)
我想念什么?
您的问题源于不使用星号取消对迭代器的引用。您将在此处的示例中注意到第28章迭代器|。每个人的Rcpp,当使用迭代器在向量的元素上循环时,为了获取特定迭代器指向的值,我们必须用星号取消引用迭代器。因此,您的代码如下所示:
for(IntegerVector::iterator i = group_index.begin(); i != group_index.end(); ++i) {
NumericVector temp = Rcpp::rbinom(m(*i) , 1, theta_sample(j, *i));
parm(0,*i) = prior_parm(0,*i) + a(*i) + sum(temp);
parm(1,*i) = prior_parm(1,*i) + b(*i) + m(*i) - sum(temp);
sample(_,*i) = Rcpp::rbeta(n_mcmc, parm(0,*i), parm(1,*i));
}
这会很快变得混乱。对我而言,最好使用范围基本循环:
for(auto i: group_index) {
NumericVector temp = Rcpp::rbinom(m(i) , 1, theta_sample(j, i));
parm(0,i) = prior_parm(0,i) + a(i) + sum(temp);
parm(1,i) = prior_parm(1,i) + b(i) + m(i) - sum(temp);
sample(_,i) = Rcpp::rbeta(n_mcmc, parm(0,i), parm(1,i));
}
有一个警告,范围基于循环。也就是说,它只是因为可用的C++11
,但是如果你使用的最新版本的不应该是一个令人关切的是大R
这使得C++11
默认。另外,您始终可以放入// [[Rcpp::plugins(cpp11)]]
源文件(请参阅将C ++ 11与Rcpp一起使用的第一步)。
另外,您正在使用(
检查向量边界的访问器,因此效率不如[
。后者更危险,并给开发人员带来责任,以确保您不访问自己不拥有的内存,但是在实践中,保证这一承诺并不困难。参见第8章第2节。rcpp给大家。
本文收集自互联网,转载请注明来源。
如有侵权,请联系 [email protected] 删除。
我来说两句