如何将 Armadillo 矩阵与从 qnorm() 获得的 NumericVector 相乘?
How can I multiply an Armadillo matrix with a NumericVector obtained from qnorm()?
当我尝试在 Rcpp 中使用 operator* 将 arma::mat 和 NumericVector 相乘时,出现以下错误:
no match for operator*
这是我要乘以的示例:
NumericVector exampleNumericVector(L-1);
arma::mat exampleMatrix(L-1,L-1);
exampleMatrix*(Rcpp::qnorm(exampleNumericVector,0,1,true,false));
我已经尝试查看 excellent Armadillo 文档,我能找到的最接近解决方案的是 prod() 函数,但这并没有解决我的问题。我也尝试过重载运算符*。然而,这也没有成功。
更新:可重现的例子
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <RcppArmadilloExtensions/sample.h>
using namespace Rcpp;
using namespace arma;
// [[Rcpp::export]]
arma::vec foo(Rcpp::NumericVector x) {
// [[Rcpp::depends(RcppArmadillo)]]
Rcpp::NumericVector y = Rcpp::qnorm(x, 0, 1, true, false);
int n = y.size();
arma::mat M(n,n);
M.fill(1.0); // nonsense content, we don't care
arma::vec v(y); // instantiate Arma vec from Rcpp vec
arma::vec res = M * v;
return res;
}
// [[Rcpp::export]]
Rcpp::List foo(const int n, const int M, const int L, const double mu){
// initialize variables
arma::mat P = arma::zeros(M,L);
arma::mat SigInvTmp = arma::zeros(L-1,L-1);
double muTmp = 0;
double s2Tmp = 0;
arma::mat Sig = arma::zeros(L,L);
arma::mat Sigee = arma::zeros(L-1,L-1);
arma::mat Sigie = arma::zeros(1,L-1);
arma::mat Sigei = arma::zeros(L-1,1);
Rcpp::NumericVector Pie(L-1);
for(int i = 0; i < M; i++){
arma::uvec iIdx(1); // vector of index to take for i
iIdx(0) = i;
for(int l = 0; l < L; l++){
srand(time(0));
// We need to get the vector of indices excluding l
arma::uvec idx(L-1); // vector of indices to subset by
arma::uword ii = 0; // the integer we'll add at each elem of idx
for (arma::uword j = 0; j < (L-1); ++j) {
if (ii == l) { // if ii equals l, we need to skip l
ii += 1;
}
idx[j] = ii; // then we store ii for this elem
ii += 1; // and increment ii
}
// Single index for l
arma::uvec lIdx(1); // vector of index to take for l
lIdx(0) = l;
Sigee = Sig.submat(idx,idx);
Sigie = Sig.submat(lIdx,idx);
Sigei = Sig.submat(idx,lIdx);
Pie = P.submat(iIdx,idx);
SigInvTmp = inv(Sigee);
arma::vec qnormVec = foo(Pie);
muTmp = mu + as_scalar(Sigie*SigInvTmp*(qnormVec-mu));
s2Tmp = sqrt(as_scalar(Sig(l,l)-Sigie*SigInvTmp*Sigei));
}
}
return Rcpp::List::create(Rcpp::Named("mu") = muTmp, Rcpp::Named("s2") = s2Tmp);
}
您的示例既不完整,也不简单,也不可重现。
所以为了批准这个,一个完整的例子。你输入的是 p 值的向量以获得分位数值,它创建一个合适维度的(无意义的)矩阵,从 Rcpp 向量创建一个 arma 向量(关键:它们都转换 to/from SEXP
)和乘以。结果返回。
代码
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::vec foo(Rcpp::NumericVector x) {
Rcpp::NumericVector y = Rcpp::qnorm(x, 0, 1, true, false);
int n = y.size();
arma::mat M(n,n);
M.fill(1.0); // nonsense content, we don't care
arma::vec v(y); // instantiate Arma vec from Rcpp vec
arma::vec res = M * v;
return res;
}
/*** R
foo(c(0.95, 0.975))
*/
演示
R> sourceCpp("~/git/Whosebug/63641766/answer.cpp")
R> foo(c(0.95, 0.975))
[,1]
[1,] 3.60482
[2,] 3.60482
R>
当我尝试在 Rcpp 中使用 operator* 将 arma::mat 和 NumericVector 相乘时,出现以下错误:
no match for operator*
这是我要乘以的示例:
NumericVector exampleNumericVector(L-1);
arma::mat exampleMatrix(L-1,L-1);
exampleMatrix*(Rcpp::qnorm(exampleNumericVector,0,1,true,false));
我已经尝试查看 excellent Armadillo 文档,我能找到的最接近解决方案的是 prod() 函数,但这并没有解决我的问题。我也尝试过重载运算符*。然而,这也没有成功。
更新:可重现的例子
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <RcppArmadilloExtensions/sample.h>
using namespace Rcpp;
using namespace arma;
// [[Rcpp::export]]
arma::vec foo(Rcpp::NumericVector x) {
// [[Rcpp::depends(RcppArmadillo)]]
Rcpp::NumericVector y = Rcpp::qnorm(x, 0, 1, true, false);
int n = y.size();
arma::mat M(n,n);
M.fill(1.0); // nonsense content, we don't care
arma::vec v(y); // instantiate Arma vec from Rcpp vec
arma::vec res = M * v;
return res;
}
// [[Rcpp::export]]
Rcpp::List foo(const int n, const int M, const int L, const double mu){
// initialize variables
arma::mat P = arma::zeros(M,L);
arma::mat SigInvTmp = arma::zeros(L-1,L-1);
double muTmp = 0;
double s2Tmp = 0;
arma::mat Sig = arma::zeros(L,L);
arma::mat Sigee = arma::zeros(L-1,L-1);
arma::mat Sigie = arma::zeros(1,L-1);
arma::mat Sigei = arma::zeros(L-1,1);
Rcpp::NumericVector Pie(L-1);
for(int i = 0; i < M; i++){
arma::uvec iIdx(1); // vector of index to take for i
iIdx(0) = i;
for(int l = 0; l < L; l++){
srand(time(0));
// We need to get the vector of indices excluding l
arma::uvec idx(L-1); // vector of indices to subset by
arma::uword ii = 0; // the integer we'll add at each elem of idx
for (arma::uword j = 0; j < (L-1); ++j) {
if (ii == l) { // if ii equals l, we need to skip l
ii += 1;
}
idx[j] = ii; // then we store ii for this elem
ii += 1; // and increment ii
}
// Single index for l
arma::uvec lIdx(1); // vector of index to take for l
lIdx(0) = l;
Sigee = Sig.submat(idx,idx);
Sigie = Sig.submat(lIdx,idx);
Sigei = Sig.submat(idx,lIdx);
Pie = P.submat(iIdx,idx);
SigInvTmp = inv(Sigee);
arma::vec qnormVec = foo(Pie);
muTmp = mu + as_scalar(Sigie*SigInvTmp*(qnormVec-mu));
s2Tmp = sqrt(as_scalar(Sig(l,l)-Sigie*SigInvTmp*Sigei));
}
}
return Rcpp::List::create(Rcpp::Named("mu") = muTmp, Rcpp::Named("s2") = s2Tmp);
}
您的示例既不完整,也不简单,也不可重现。
所以为了批准这个,一个完整的例子。你输入的是 p 值的向量以获得分位数值,它创建一个合适维度的(无意义的)矩阵,从 Rcpp 向量创建一个 arma 向量(关键:它们都转换 to/from SEXP
)和乘以。结果返回。
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::vec foo(Rcpp::NumericVector x) {
Rcpp::NumericVector y = Rcpp::qnorm(x, 0, 1, true, false);
int n = y.size();
arma::mat M(n,n);
M.fill(1.0); // nonsense content, we don't care
arma::vec v(y); // instantiate Arma vec from Rcpp vec
arma::vec res = M * v;
return res;
}
/*** R
foo(c(0.95, 0.975))
*/
演示
R> sourceCpp("~/git/Whosebug/63641766/answer.cpp")
R> foo(c(0.95, 0.975))
[,1]
[1,] 3.60482
[2,] 3.60482
R>