The rate/scale parameter in rexp

[xiangjiexue]

Hi,

The parameter of rexp() is the rate in TMB. I went to check the source code of the implementation, this function calls the C implementation in R, which is essentially the function located in /src/nmath/rexp.c. I have checked the source code of this file (at least R 4.0.4), the function actually takes the scale as the parameter.

Here is a minimal example estimating the rate parameter of the exponential distribution and performs a simulation on the estimated parameter.

C++ side:

#include <TMB.hpp>

template<class Type>
Type objective_function<Type>::operator() ()
{
  using namespace density; 
  using namespace Eigen; 

  DATA_VECTOR(y); 
  PARAMETER(log_lambda);
  Type lambda = exp(log_lambda);

  Type nll = -sum(dexp(y, lambda, true));

  SIMULATE {
  	int L = y.size();
  	for (int i = 0; i < L; ++i){
  		y[i] = rexp(lambda);
  	}
  	REPORT(y);
  }

  return nll;
}

and the R side

library(TMB)

lambda = 10
y = rexp(1000, lambda)

compile("rexptest.cpp")
dyn.load("rexptest.so")

data = list(y = y)
param = list(log_lambda = 0)
obj = MakeADFun(data, param, silent = TRUE, hessian = TRUE)
opt <- nlminb(obj$par,obj$fn,obj$gr) // This give 2.26 for log of lambda.

sim <- replicate(100, {
  simdata <- obj$simulate(complete=TRUE)
  obj2 <- MakeADFun(simdata, param, silent=TRUE)
  nlminb(obj2$par, obj2$fn, obj2$gr)$par
})
mean(sim) // This gives around -2.26, but log(10) = 2.30.

dyn.unload("rexptest.so")

Thanks,
Xiangjie