I'm a newbie with this, so sorry for the basic question. I have managed to get Rcpp working but need a function in Armadillo, and I am just not getting it right. I know the code is bad, like I said I am a newbie in C++, but I would be very grateful if someone could help me get this example going. My code from the .cpp file in Rstudio is as follows. I am running latest R, and Rstudio. Thanks.
#include <iostream>
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;
/*** R
require(Rcpp)
require(RcppArmadillo)
#Population of size N
N = 350
#Set size
m = 3
#Number of cycles
C = 20
#Final number of analysed units
n = C * m
#make the r vector of ranks
r <- rep(1:m, length.out = n)
#initialise alpha matrix nxn
alpha <- matrix(0, ncol = N, nrow = n)
*/
// [[Rcpp::export]]
unsigned nChoosek( unsigned n, unsigned k )
{
if (k > n) return 0;
if (k * 2 > n) k = n-k;
if (k == 0) return 1;
int result = n;
for( int i = 2; i <= k; ++i ) {
result *= (n-i+1);
result /= i;
}
return result;
}
// [[Rcpp::export]]
NumericMatrix alphaC(NumericVector r, int N, int m, int n, int C,
NumericMatrix alpha){
int S = r.size();
int a = 0;
int b = 0;
double c = 0;
for(int i = 0; i < N; ++i){
for (int j = 0; j < S; ++j){
a = nChoosek(i, (r[j] - 1));
b = nChoosek((N - (i + 1)),( m - r[j]));
c = nChoosek(N, m);
alpha(j,i) = a*b/c;
//std::cout << "a*b/c: " << a*b/c << std::endl;
}
}
return(alpha);
}
/*** R
#Create first order probabilities
pi <- 1 - apply((1 - alphaC(r,N,m,n,C,alpha)), 2, prod)
plot(pi)
*/
/*** R
#Second order inclusion properties.
pii <- matrix(0, ncol = N, nrow = N)
*/
// [[Rcpp::export]]
int times(NumericVector vec){
int sum = 1;
sum = arma::prod(vec);
return sum;
}
/*** R
vec = 1:5
times(vec)
*/
Your example code contains many things that are not relevant here. Please try to provide a minimal example next time, together with the error message you got. Anyway, you will need an Armadillo data structure to work with arma::prod
. You can either create this directly as shown by @mt1022 in the comments, or have it done implicitly by defining a function with suitable arguments:
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
double times(arma::vec vec){
double result = arma::prod(vec);
return result;
}
/*** R
vec <- 1:5 + 0.5
times(vec)
*/
Note that I have changed the return type to double, since the product of a numeric
vector is in general not representable as a int
.