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

Question

我在 r that is running quite slowly, and so I'm hoping to rewrite the function in 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
    })

此代码应 return 嵌套 list of matrices. The list should be of length numstudies and each element of it should be a matrix 维度 n[v] 行乘 ql 列。这些矩阵的每一行都应该从多元正态分布中抽取，但是每行的抽取将由不同的均值向量 finalmu_b2_l 和不同的协方差矩阵 finalcov_b2_l 控制。 finalmu_b2_l和finalcov_b2_l都是嵌套列表，顶层长度numstudies，每个元素都是一个长度n[v].

的列表

到目前为止，我已经尝试编写一段 rcpp 代码来绘制所需的矩阵之一，但想尽快将其扩展为 return 矩阵列表可能的。代码如下：

#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 文件，使用函数 sourceCpp），我都会收到错误消息：

invalid initialization of reference of type 'const::mat&'

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

关于 -

的任何建议、指示或指导

为什么这段代码没有编译，我应该在 Rcpp 中做什么，and/or
如何将其扩展为 return 矩阵列表而不是单个矩阵将不胜感激。

编辑

原始 r 代码包含嵌套在 lapply 函数中的 mapply 函数。我想要一个列表长度 numstudies 的输出，其中的每个元素都是一个矩阵。此列表中每个矩阵的每一行都是来自不同多元正态分布的 reliation（因此在每个列表元素中，对于矩阵的每一行，都有一个唯一的均值向量和一个唯一的协方差矩阵控制我想要绘制的多元正态分布从）。我写成一个嵌套在 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;
}

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

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() 中，正如您所说，它们应该是列表的列表。）

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

Rcpp function with List arguments returning a matrix

c++

performance

r

rcpp

rcpparmadillo