带有返回矩阵的List参数的Rcpp函数

Question

我在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;
}

每当我试图编译代码时(通过使用函数sourceCpp获取代码所在的.cpp文件(，我都会得到一个错误：

类型为"const:：mat&"的引用初始化无效

，我认为这是创建draws矩阵或填充矩阵的问题？

关于-的任何建议、指示或指导

1. 为什么不编译此代码，我应该在Rcpp中做什么，和/或
2. 如何将其扩展为返回矩阵列表而不仅仅是单个矩阵将不胜感激

编辑

原始r代码包含一个嵌套在lapply函数中的mapply函数。我想要一个输出，它是一个列表长度的numstudies，其中的每个元素都是一个矩阵。该列表中的每个矩阵的每一行都是来自不同的多元正态分布的重化(因此，在每个列表元素中，对于矩阵的每行，都有一个唯一的均值向量和一个唯一协方差矩阵来控制我想要从中提取的多元正模(。我编写了一个嵌套在lapply中的mapply，以自动给出我想要的输出格式，同时允许为每个矩阵行从不同的分布中绘制。

EDIT2

将迭代更改为从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;
}

EDIT3

该代码当前正在编译，但未绘制示例。相反，我得到以下错误消息：

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;
}

当然，一旦它起作用，我仍然需要扩展它来绘制矩阵列表，而不是单个矩阵

Answer 1

这里有一个基本的设置，应该可以做你想要的：

#include <RcppDist.h>
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppDist,RcppArmadillo)]]
// [[Rcpp::export]]
arma::mat draw_b(const int ql,
const Rcpp::List& mu,
const Rcpp::List& cov) {
int n = mu.size();
arma::mat draws = arma::zeros(n, ql);
for ( arma::uword iter = 0; iter < n; iter++ ) {
draws.row(iter) = rmvnorm(1, mu[iter], cov[iter]);
}
return draws;
}
// [[Rcpp::export]]
Rcpp::List get_list_of_draws(Rcpp::List mu, Rcpp::List Sigma, int ql) {
int numstudies = mu.size();
Rcpp::List res(numstudies);
for ( int iter = 0; iter < numstudies; ++iter ) {
Rcpp::List mu_temp = mu[iter];
Rcpp::List cov_temp = Sigma[iter];
res[iter] = draw_b(ql, mu_temp, cov_temp);
}
return res;
}

它似乎如预期的那样工作：

res <- get_list_of_draws(mu_example, cov_example, ql)
str(res)
# List of 1
#  $ : num [1:10, 1:2] -0.0156 -0.4717 -0.8134 0.5489 0.1215 ...

(不过请注意，当我设置mu_example和cov_example时，我将它们封装在list()中，正如您所说，它们应该是列表列表。(