使用 RcppArmadillo::sample 增量采样时的奇怪行为
Strange behavior when incrementally sampling using RcppArmadillo::sample
我正在尝试使用 Rcpp 使用 polya urn 方案实现一些绘图。基本上,我有一个正在绘制的矩阵,以及一个权重与概率成正比的第二矩阵。每次绘制后,我需要增加我绘制的单元格的权重。
我 运行 遇到了一些索引错误,这导致我更广泛地检查抽样,我发现我的权重矩阵被 RcppArmadillo::sample 修改了。两个问题 (1) 这种行为是我应该预料到的还是我应该在某处报告的错误? (2) 关于当前解决方法的任何想法?这是一个可重现的例子:
#include <RcppArmadilloExtensions/sample.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp ;
// [[Rcpp::export]]
void sampler(int N, int inc, NumericMatrix& weight_matrix, int reps) {
IntegerVector wm_tmp = seq_along(weight_matrix);
Rcout << "Initial weight_matrix:\n" << weight_matrix << "\n";
int x_ind;
for(int i = 0; i < reps; ++i) {
x_ind = RcppArmadillo::sample(wm_tmp, 1, true, weight_matrix)(0) - 1;
Rcout << "Weight matrix after sample: (rep = " << i << ")\n" << weight_matrix << "\n";
Rcout << "x_ind: " << x_ind << "\n";
// get indices
weight_matrix[x_ind] = weight_matrix[x_ind] + inc;
Rcout << "Add increment of " << inc << " to weight_matrix:\n" << weight_matrix << "\n";
}
}
//
// // [[Rcpp::export]]
// IntegerVector seq_cpp(IntegerMatrix x) {
// IntegerVector tmp = seq_along(x);
// IntegerVector ret = RcppArmadillo::sample(tmp, 2, true);
// return ret;
// }
/*** R
weight_matrix <- matrix(1, 5, 2)
sampler(5, 1, weight_matrix, 3)
weight_matrix <- matrix(1, 5, 2)
sampler(5, 0, weight_matrix, 3)
*/
谢谢!
这是已知且有记录的行为。
你可以
i) 使用 Rcpp::clone()
创建您的 SEXP
的不同副本(即 NumericMatrix
)。
ii) 使用 Armadillo 矩阵代替并传递为 const arma::mat & m
.
存在与 R 组织其数据结构的方式有关的体系结构原因,这意味着我们无法为您提供快速访问(无副本!)并防止写入。
我正在尝试使用 Rcpp 使用 polya urn 方案实现一些绘图。基本上,我有一个正在绘制的矩阵,以及一个权重与概率成正比的第二矩阵。每次绘制后,我需要增加我绘制的单元格的权重。
我 运行 遇到了一些索引错误,这导致我更广泛地检查抽样,我发现我的权重矩阵被 RcppArmadillo::sample 修改了。两个问题 (1) 这种行为是我应该预料到的还是我应该在某处报告的错误? (2) 关于当前解决方法的任何想法?这是一个可重现的例子:
#include <RcppArmadilloExtensions/sample.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp ;
// [[Rcpp::export]]
void sampler(int N, int inc, NumericMatrix& weight_matrix, int reps) {
IntegerVector wm_tmp = seq_along(weight_matrix);
Rcout << "Initial weight_matrix:\n" << weight_matrix << "\n";
int x_ind;
for(int i = 0; i < reps; ++i) {
x_ind = RcppArmadillo::sample(wm_tmp, 1, true, weight_matrix)(0) - 1;
Rcout << "Weight matrix after sample: (rep = " << i << ")\n" << weight_matrix << "\n";
Rcout << "x_ind: " << x_ind << "\n";
// get indices
weight_matrix[x_ind] = weight_matrix[x_ind] + inc;
Rcout << "Add increment of " << inc << " to weight_matrix:\n" << weight_matrix << "\n";
}
}
//
// // [[Rcpp::export]]
// IntegerVector seq_cpp(IntegerMatrix x) {
// IntegerVector tmp = seq_along(x);
// IntegerVector ret = RcppArmadillo::sample(tmp, 2, true);
// return ret;
// }
/*** R
weight_matrix <- matrix(1, 5, 2)
sampler(5, 1, weight_matrix, 3)
weight_matrix <- matrix(1, 5, 2)
sampler(5, 0, weight_matrix, 3)
*/
谢谢!
这是已知且有记录的行为。
你可以
i) 使用 Rcpp::clone()
创建您的 SEXP
的不同副本(即 NumericMatrix
)。
ii) 使用 Armadillo 矩阵代替并传递为 const arma::mat & m
.
存在与 R 组织其数据结构的方式有关的体系结构原因,这意味着我们无法为您提供快速访问(无副本!)并防止写入。