Использование rmultinom с Rcpp

Я хотел бы использовать функцию R rmultinom в коде c ++ для использования с Rcpp. Я получаю ошибку о недостаточном количестве аргументов — я не знаю, какими должны быть эти аргументы, поскольку они не соответствуют аргументам, используемым с функцией в R. Мне также не повезло с использованием «:: Rf_foo» синтаксис для доступа к функции R из кода Rcpp.

Вот упрощенная версия моего кода ниже (да, я пишу сэмплер Гиббса).

#include <Rcpp.h>
using namespace Rcpp;

// C++ implementation of the R which() function.
int whichC(NumericVector x, double val) {
int ind = -1;
int n = x.size();
for (int i = 0; i < n; ++i) {
if (x[i] == val) {
if (ind == -1) {
ind = i;
} else {
throw std::invalid_argument( "value appears multiple times." );
}
} // end if
} // end for
if (ind != -1) {
return ind;
} else {
throw std::invalid_argument( "value doesn't appear here!" );
return -1;
}
}

// [[Rcpp::export]]
int multSample(double p1, double p2, double p3) {
NumericVector params(3);
params(0) = p1;
params(1) = p2;
params(2) = p3;

// HERE'S THE PROBLEM.
RObject sampled = rmultinom(1, 1, params);
int out = whichC(as<NumericVector>(sampled), 1);
return out;
}

Я новичок в C ++, так что я понимаю, что большая часть этого кода, вероятно, Nooby и неэффективно. Я открыт для предложений о том, как улучшить мой код на C ++, но мой приоритет — услышать о бизнесе rmultinom. Спасибо!

Кстати, я прошу прощения за сходство с эта тема, но

  1. Ответ не работал для моих целей
  2. Разницы может быть достаточно, чтобы оправдать другой вопрос (вы так думаете?)
  3. Вопрос был опубликован и ответил год назад.

3

Решение

Ниже приведен ответ пользователя user95215 с поправкой на его компиляцию и еще одна версия в стиле Rcpp:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
IntegerVector oneMultinomC(NumericVector probs) {
int k = probs.size();
SEXP ans;
PROTECT(ans = Rf_allocVector(INTSXP, k));
probs = Rf_coerceVector(probs, REALSXP);
rmultinom(1, REAL(probs), k, &INTEGER(ans)[0]);
UNPROTECT(1);
return(ans);
}

// [[Rcpp::export]]
IntegerVector oneMultinomCalt(NumericVector probs) {
int k = probs.size();
IntegerVector ans(k);
rmultinom(1, probs.begin(), k, ans.begin());
return(ans);
}
3

Другие решения

Если я пытаюсь скомпилировать ваш код, я получаю ошибку компилятора:

> Rcpp::sourceCpp('~/scratch/multSample.cpp')
multSample.cpp:33:21: error: no matching function for call to 'rmultinom'
RObject sampled = rmultinom(1, 1, params);
^~~~~~~~~
/Library/Frameworks/R.framework/Resources/include/Rmath.h:449:6: note: candidate function not viable: requires 4 arguments, but 3 were provided
void    rmultinom(int, double*, int, int*);
^
1 error generated.

Как следует из вышесказанного, вы не указали аргументы правильно. Обратите внимание, что rmultinom интерфейс немного неловкий по сравнению с другими функциями: он заполняет память, на которую указывает *rnвместо того, чтобы возвращать новый объект (с его собственной недавно выделенной памятью).

Если вы посмотрите в источники R вы увидите интерфейс, и вы можете увидеть пример его использования Вот а также (на самом деле, stats делает функцию-обертку, которая еще проверяет аргументы и еще много чего). Но обратите внимание, как это используется здесь:

rmultinom(size, REAL(prob), k, &INTEGER(ans)[ik]);

Другими словами, это заполняет INTSXP называется ans передав этот адрес этой памяти на rmultinom функция.

Так что если вы хотите использовать эту функцию из Rcppвам придется сделать что-то подобное — но, возможно, это заслуживает аналогичной обработки векторизации сахара, чтобы избежать уродства этого интерфейса.

Вы можете попробовать сделать что-то вроде:

IntegerMatrix sampled(nrow, ncol);
rmultinom(1, 1, params, sampled.begin());

Или что-то в этом роде.

0

Пример и ссылки, предоставленные Кевином, позволили мне найти ответ, который сработал. Был какой-то спор типов. Я написал функцию, которая позволяет вам выбирать один вектор из полиномиального распределения. Код ниже.

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector oneMultinomC(NumericVector probs) {
int k = probs.size();
SEXP ans;
PROTECT(ans = RF_allocVector(INTSXP, k));
probs = RF_coerceVector(probs, REALSXP);
rmultinom(1, REAL(probs), k, &INTEGER(ans)[0]);
UNPROTECT(1);
return ans;
}

Я не понимаю половину того, что здесь происходит. В частности, я не понимаю четвертого аргумента к «рмультином». Я знаю, что это указатель на место в памяти, где хранится вывод, но я не понимаю бит «[0]». Тем не менее, это работает. Образцы Гиббса далеко, мальчики и девочки.

0