Как я могу умножить матрицу броненосца на числовой вектор, полученный из qnorm()?

Когда я пытаюсь умножить arma::mat и NumericVector в Rcpp с помощью оператора*, я получаю следующую ошибку:

no match for operator*

Вот пример того, что я пытаюсь умножить:

NumericVector exampleNumericVector(L-1);
arma::mat exampleMatrix(L-1,L-1);
exampleMatrix*(Rcpp::qnorm(exampleNumericVector,0,1,true,false));

Я попытался просмотреть отличную документацию 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);
}

person Carl Yocto    schedule 28.08.2020    source источник
comment
Вероятно, проще всего преобразовать матрицу в NumericMatrix или вектор в матрицу броненосца.   -  person Oliver    schedule 29.08.2020
comment
@Oliver Я рассматривал оба этих варианта, отдавая предпочтение последнему. Однако я не могу найти способ преобразовать NumericVector в матрицу Armadillo.   -  person Carl Yocto    schedule 29.08.2020
comment
@Oliver Приведение типов одинаково для перехода от numericVector к arma::mat? не могу найти функцию для этого   -  person Carl Yocto    schedule 30.08.2020
comment
Привет @RonSnow. Ответ Дирка - это прямой метод вашего конкретного решения. В С++ приведение типов выполняется с использованием (new-type)variable возможного с указателем ссылки. Для этого требуется, чтобы был определен метод, который может привести метод от одного к другому. Его ответ гораздо более общий и удобный. В то время как Дирк может быть очень прямолинейным в своем отношении (в основном из-за количества ответов, которые он дает, я думаю), его ответы блестящие. Для вашей конкретной проблемы я бы переместил // [[Rcpp::depends(RcppArmadillo)]] к фактической функции, хотя я не уверен, что это ответ.   -  person Oliver    schedule 30.08.2020
comment
@Oliver Спасибо за внимание. Я реализовал функцию Дирка, но по-прежнему получаю сообщение об ошибке статического утверждения: невозможно преобразовать тип в SEXP. Я думаю, что это результат моей строки кода muTmp = mu + as_scalar(Sigie*SigInvTmp*(qnormVec-mu)); Как вы думаете?   -  person Carl Yocto    schedule 30.08.2020


Ответы (1)


Ваш пример не является ни полным, ни минимальным, ни воспроизводимым.

Итак, в попытке утвердить это, полный пример. Вы вводите вектор значений p для получения значений квантилей, он создает (бессмысленную) матрицу подходящей размерности, создает вектор arma из вектора Rcpp (ключ: все они преобразуются в/из SEXP) и умножаются. Результат возвращается.

Code
#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))
*/
Demo
R> sourceCpp("~/git/stackoverflow/63641766/answer.cpp")

R> foo(c(0.95, 0.975))
        [,1]
[1,] 3.60482
[2,] 3.60482
R> 
person Dirk Eddelbuettel    schedule 29.08.2020
comment
Спасибо за ваш пример. Я попытался применить ваш метод к своей проблеме, и я получаю следующую ошибку: статическое утверждение не удалось: невозможно преобразовать тип в SEXP. - person Carl Yocto; 29.08.2020
comment
Пожалуйста, прочитайте stackoverflow.com/help/minimal-reproducible-example -- наше устройство для чтения мыслей находится в магазине, поэтому мы не могу сказать, какой код вы написали. - person Dirk Eddelbuettel; 29.08.2020
comment
Мои извинения - у вас есть хорошая мысль. Я обновил сообщение, включив в него минимальный воспроизводимый пример. - person Carl Yocto; 29.08.2020