Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Poisson draw in Rcpp and R different results

Tags:

r

rcpp

I face the following contradiction when I use the same code in R and in Rcpp

In R I run the following code

t = 0
for(i in 1:50){
 t = t + rpois(1, 0.5)
}
t
[1] 28

and I take back a value t which is nonnegative. Now I type the exact same commands in Rcpp

#include <Rcpp.h>
#include<Rmath.h>
using namespace Rcpp;

// [[Rcpp::export]]
int Pois(int l){ 
  int t=0;
  for(int i=0; i<50;++i){
    t+=R::rpois(l);
  }
  return t;
}

and when I call the function in R

Pois(0.5)
[1] 0

which is wrong since in R it was different of zero

What is going wrong?

like image 707
Jonathan1234 Avatar asked Jan 31 '26 02:01

Jonathan1234


2 Answers

You should use double l rather than int l, e.g.,

int Pois(double l){ 
  int t=0;
  for(int i=0; i<50;++i){
    t+=R::rpois(l);
  }
  return t;
}

otherwise (int) 0.5 gives you 0.

like image 136
ThomasIsCoding Avatar answered Feb 02 '26 17:02

ThomasIsCoding


@ThomasIsCoding already showed you the main issue. But recall that beside R::rpois() we also have the vectorised Rcpp::rpois(). And, as usual, given the same seed it gives the same draws as R:

> set.seed(123)
> rpois(10, 0.5)
 [1] 0 1 0 1 2 0 0 1 0 0
> Rcpp::cppFunction("NumericVector myrp(int n, double l) { return Rcpp::rpois(n, l); }")
> set.seed(123)
> myrp(10, 0.5)
 [1] 0 1 0 1 2 0 0 1 0 0
> 
like image 35
Dirk Eddelbuettel Avatar answered Feb 02 '26 19:02

Dirk Eddelbuettel



Donate For Us

If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!