我在r中有一段代码运行得很慢,所以我希望重写rcpp中的函数,但是,这是我第一次尝试使用 Rcpp,我无法编译或运行代码。
在 R 中,我试图重写的部分是
drawsamples<-lapply(1:numstudies,function(v){
temp<-t(mapply(rmvn,n=1,mu=finalmu_b2_l[[v]],sigma=finalcov_b2_l[[v]]))
rownames(temp)<-id[[v]]
temp
})
此代码应返回一个嵌套的矩阵列表。该列表应该是长度numstudies
,并且它的每个元素应该是一个按列排列的维度矩阵。这些矩阵的每一行都应该从多元正态分布中抽取,但是每一行的抽取将由不同的均值向量和不同的协方差矩阵控制。两者和都是嵌套列表,顶级长度,其中每个元素都是长度列表。n[v]
ql
finalmu_b2_l
finalcov_b2_l
finalmu_b2_l
finalcov_b2_l
numstudies
n[v]
到目前为止,我已经尝试编写一段rcpp代码,它只绘制一个所需的矩阵,但想扩展它以尽快返回一个矩阵列表。代码如下:
#include <RcppDist.h>
#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::depends(RcppDist,RcppArmadillo)]]
// [[Rcpp::export]]
arma::mat draw_b(const int n,
const int ql,
const Rcpp::List mu,
const Rcpp::List cov) {
arma::mat draws = arma::zeros(n,ql);
for (int iter = 1; iter <= n; iter++) {
draws.row(iter) = rmvnorm(1, mu[iter], cov[iter]);
}
return draws;
}
每当我尝试编译代码时(通过使用 .cpp 文件获取代码所在的位置,使用 function sourceCpp
),我都会收到错误消息:
'const::mat&' 类型的引用初始化无效
,我将其理解为创建draws
矩阵或填充矩阵的问题?
任何建议、指示或指导 -
- 为什么这段代码没有编译,我应该在 Rcpp 中做什么,和/或
- 如何扩展它以返回矩阵列表而不仅仅是单个矩阵将不胜感激。
编辑
原始 r 代码包含一个嵌套在 lapply 函数中的 mapply 函数。我想要一个输出是一个列表长度 numstudies,其中的每个元素都是一个矩阵。该列表中每个矩阵的每一行都是来自不同的多元正态分布的关联(因此在每个列表元素中,对于矩阵的每一行,都有一个唯一的均值向量和一个唯一的协方差矩阵来控制我想要绘制的多元正态分布从)。我编写了一个嵌套在 lapply 中的 mapply,以自动给出我想要的输出格式,同时允许从每个矩阵行的不同分布中进行绘制。
编辑2
将迭代更改为从 0 而不是 1 运行后,将编译以下代码:
#include <RcppDist.h>
#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::depends(RcppDist,RcppArmadillo)]]
//' @keywords internal
// [[Rcpp::export]]
arma::mat draw_b(const int & n,
const int & ql,
const Rcpp::List & mu,
const Rcpp::List & cov) {
arma::mat draws = arma::zeros(n,ql);
for (int iter = 0; iter <= n; iter++ ) {
draws.row(iter) = rmvnorm(1, mu_temp, cov_temp);
}
return draws;
}
编辑3
该代码当前编译,但不绘制样本。相反,我收到以下错误消息:
error: Mat::init(): requested size is not compatible with column vector layout
Error in draw_b(n = n, ql = ql, mu = mu_example, cov = cov_example) :
Mat::init(): requested size is not compatible with column vector layout
我已经准备了一个我希望这个基本函数做什么的例子(它只是从不同的多元正态分布中采样一个实现矩阵。
数据:
list(n = 10, ql = 2, mu_example = list(c(0.342909965003207, -0.788070875792469
), c(-0.00499810116271365, 0.0114865660452949), c(-0.149753928200309,
0.344162379034614), c(0.335176829763227, -0.770298692761465),
c(0.254206123984596, -0.584212951520601), c(0.379893097582703,
-0.873064992779416), c(0.137231089484867, -0.315382566602526
), c(0.405123380985852, -0.931048876501857), c(-0.00505917031396947,
0.0116269143128456), c(-0.0743318653279181, 0.170828451158346
)), cov_example = list(structure(c(0.199912910315971, -0.459437048770092,
-0.459437048770092, 4.49223135519527), .Dim = c(2L, 2L)), structure(c(0.199912910315971,
-0.459437048770092, -0.459437048770092, 4.49223135519527), .Dim = c(2L,
2L)), structure(c(0.199912910315971, -0.459437048770092, -0.459437048770092,
4.49223135519527), .Dim = c(2L, 2L)), structure(c(0.199912910315971,
-0.459437048770092, -0.459437048770092, 4.49223135519527), .Dim = c(2L,
2L)), structure(c(0.199912910315971, -0.459437048770092, -0.459437048770092,
4.49223135519527), .Dim = c(2L, 2L)), structure(c(0.199912910315971,
-0.459437048770092, -0.459437048770092, 4.49223135519527), .Dim = c(2L,
2L)), structure(c(0.199912910315971, -0.459437048770092, -0.459437048770092,
4.49223135519527), .Dim = c(2L, 2L)), structure(c(0.199912910315971,
-0.459437048770092, -0.459437048770092, 4.49223135519527), .Dim = c(2L,
2L)), structure(c(0.199912910315971, -0.459437048770092, -0.459437048770092,
4.49223135519527), .Dim = c(2L, 2L)), structure(c(0.199912910315971,
-0.459437048770092, -0.459437048770092, 4.49223135519527), .Dim = c(2L,
2L))))
R代码
drawsampletest<-draw_b(n=n,
ql=ql,
mu=mu_example,
cov=cov_example)
rcpp 代码
#include <RcppDist.h>
#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::depends(RcppDist,RcppArmadillo)]]
//' @keywords internal
// [[Rcpp::export]]
arma::mat draw_b(const int & n,
const int & ql,
const Rcpp::List & mu,
const Rcpp::List & cov) {
arma::mat draws = arma::zeros(n,ql);
for (int iter = 0; iter <= n; iter++ ) {
arma::rowvec mu_temp = mu[iter];
arma::mat cov_temp = cov[iter];
draws.row(iter) = rmvnorm(1, mu_temp, cov_temp);
}
return draws;
}
当然,一旦它起作用了-我仍然需要扩展它以绘制矩阵列表而不是单个矩阵